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ABSTRACT 

Recently, we explored new, meshless finite-volume Lagrangian methods for hydrodynamics: the “meshless finite 
mass” (MFM) and “meshless finite volume” (MFV) methods; these capture advantages of both smoothed-particle 
hydrodynamics (SPH) and adaptive mesh-refinement (AMR) schemes. We extend these to include ideal magnetohy¬ 
drodynamics (MHD). The MHD equations are second-order consistent and conservative. We augment these with a 
divergence-cleaning scheme, which maintains V ■ B ss 0. We implement these in the code GIZMO, together with state- 
of-the-art SPH MHD. We consider a large test suite, and show that on all problems the new methods are competitive 
with AMR using constrained transport (CT) to ensure V ■ B = 0. They correctly capture the growth/structure of the 
magnetorotational instability (MRI), MHD turbulence, and launching of magnetic jets, in some cases converging more 
rapidly than state-of-the-art AMR. Compared to SPH, the MFM/MFV methods exhibit convergence at fixed neighbor 
number, sharp shock-capturing, and dramatically reduced noise, divergence errors, & diffusion. Still, “modern” SPH 
can handle most test problems, at the cost of larger kernels and “by hand” adjustment of artificial diffusion. Compared 
to non-moving meshes, the new methods exhibit enhanced “grid noise” but reduced advection errors and diffusion, 
easily include self-gravity, and feature velocity-independent errors and superior angular momentum conservation. 
They converge more slowly on some problems (smooth, slow-moving flows), but more rapidly on others (involving 
advection/rotation). In all cases, we show divergence-control beyond the Powell 8-wave approach is necessary, or all 
methods can converge to unphysical answers even at high resolution. 
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1 INTRODUCTION: THE CHALLENGE OF EXISTING 
NUMERICAL METHODS 

Magnetic fields are an essential component in astrophysical hy¬ 
drodynamics, and for many astrophysical problems can be reason¬ 
ably approximated by ideal (infinite conductivity) magnetohydro¬ 
dynamics (MHD). The MHD equations are inherently non-linear, 
however, and most problems require numerical simulations. But 
this poses unique challenges, especially for numerical methods 
which are Lagrangian (i.e. the mesh elements follow the fluid), 
rather than Eulerian (solved on a fixed grid). 

In most discretizations (although see |Kawai|2013[ >, evolving 
the MHD equations in time will lead to a violation of the “diver¬ 
gence constraint” (the requirement that V ■ B = 0). Unfortunately, 
this cannot simply be ignored to treated as a “standard” numeri¬ 
cal error term which should converge away with increasing reso¬ 
lution, because certain errors introduced by a non-zero V ■ B are 
numerically unstable: they will eventually destroy the correct so¬ 
lution (even at infinite resolution) and/or produce unphysical re¬ 
sults (e.g. negative pressures). Arguably the most elegant solution 
is the so-called “constrained transport” (CT) method of [Evans~&] 
|Hawley| ( |1988) , which maintains V ■ B to machine precision; how¬ 
ever, while there is no obvious barrier in principle to implementing 
this in meshless and unstructured mesh methods (see recent devel¬ 
opments by |Mocz et al.||2014a| l, it has thus far only been practi¬ 
cal to implement for real problems in regular, Cartesian grid (or 
adaptive-mesh refinement; AMR) codes. But for many problems in 
astrophysics, Lagrangian, mesh-free codes have other advantages: 
they minimize numerical diffusion and over-mixing, move with the 
fluid so automatically provide enhanced resolution with the mass 
(in a continuous manner, which avoids low-order errors necessar¬ 
ily introduced by AMR refinement boundaries), couple simply and 
accurately to cosmological expansion and A-body gravity codes, 
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easily handle high Mach numbers, conserve angular momentum 
and naturally handle orbiting disks without prior knowledge of the 
disk geometry, avoid “grid alignment” and carbuncle instabilities 
(where the grid imprints preferred directions on the gas), and fea¬ 
ture errors which are independent of the fluid bulk velocity (so can 
converge more rapidly when the fluid moves). 

A variety of approaches have been developed to deal with 
these errors. The simplest commonly-used method, the so-called 
“Powell 8-wave cleaning,” simply subtracts the unstable error terms 
resulting from a non-zero V ■ B from the equation of motion (Pow- 
|ell et al.|1999) . This removes the more catastrophic numerical in¬ 
stabilities, but does not solve the convergence problem - many 
studies have shown that certain types of problems, treated with 
only this method, will simply converge to the wrong solution (jTothj 
|2000| |Mignone & Tzeferacos|2010{|Mocz et al,|2QI4aj . And the 
subtraction necessarily violates momentum conservation, so one 
would ideally like the subtracted terms (the V ■ B values) to re¬ 
main as small as possible. Therefore more sophisticated “cleaning” 
schemes have been developed, the most popular of which have been 
variants of the |Dedner et al.| ( f2002]> method: this adds source terms 
which transport the divergence away (in waves) and then damp it. 
This has proven remarkably robust and stable. 

However, applications of these techniques in Lagrangian 
codes in astrophysics have remained limited. The most popular La¬ 
grangian method, smoothed-particle hydrodynamics (SPH), suffers 
from several well-known errors that make MHD uniquely challeng¬ 
ing. The SPH equations are not consistent at any order (meaning 
they contain zeroth-order errors; |Morris|[T996| |Dilts||1999{ [Read] 
|et al.||2010| l; this introduces errors which converges away very 
slowly and causes particular problems for divergence-cleaning. 
Also, naive implementations of the equations are vulnerable to 
the tensile and particle pairing instabilities. And artificial diffusion 
terms, with ad-hoc parameters, are required in SPH to deal with 
discontinuous fluid quantities. As such, many previous implemen¬ 
tations of MHD in SPH were unable to reproduce non-trivial field 
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configurations, were extremely diffusive, or were simply unable to 
numerically converge; in turn key qualitative phenomena such as 
the magneto-rotational instability (MRI) and launching of magnetic 
jets could not be treated (see |Swegle_et_alj 1995; Monaghan 2000 


B0rvee^alJ2OOljjMaron&HowesJ2OO3jJPrice^Rosswog 2006 

Rosswog & Price[|2007| |Price & Bate||2008[ |Dolag & Stasyszyn| 


2009^ 


Recently, however, a number of breakthroughs have been 
made in Lagrangian hydrodynamics, with the popularization of 
moving-mesh and mesh-free finite-volume Godunov methods. 
|Springel| {2010}; |Duffell & MacFadyen| l |2011| >; |Gaburov et al.| 
(2012) have developed moving-mesh MHD codes, which capture 


many of the advantages of both AMR and SPH, using the Ded- 


|ner et ah] ( 2002} cleaning meth od. Meanwhile, |Lanson & Vila 
1 2008a|b i; |Gaburov & Nitadori| ( [2011| >; |Hopkins| 1 2015} have de 


veloped a class of new, mesh-free finite volume methods which 
are both high-order consistent (convergent) and fully conservative. 
These are very similar to moving-mesh codes (in fact, Voronoi 
moving-meshes are technically a special case of the method). In 
|Hopkins|2015[ these are developed for hydrodynamics in the multi¬ 
method, hydrodynamics+gravity+cosnrology code GIZMO, which 
is an extension of the A-body gravity and domain decomposition 
algorithms from GADGET-3 (Springe f|2005} to i nclude a variety of 

s0I 


new hydrodynamic methodsrjln Hopkins (20151, a broad range of 


test problems are considered, and it is shown that these also capture 
most of the advantages of AMR and SPH, while avoiding many of 
their disadvantages. Particularly important, these eliminate the low- 
order errors, numerical instabilities, and artificial diffusion terms 
which have plagued SPH. Gaburov & Nitadori (2011) considered a 
range of MHD test problems and found very encouraging prelimi¬ 
nary results; they showed that these mesh-free methods could han¬ 
dle complicated non-linear problems like the MRI with accuracy 
comparable to state-of-the-art grid codes. Meanwhile, tremendous 


2001; Pricej2008 

jWadsley et al.||2008[ Cullen & Dehnen 2010 

Read & Hayfield 

2012, Saitoh & Makino 2013, Hopkins 2013 


Therefore, in this paper, we extend the nresh-free MFM and 
MFV Lagrangian hydrodynamics in GIZMO to include MHD, and 
consider a systematic survey of a wide range of test problems, and 
compare state-of-the-art grid-based (AMR) codes, MFM, MFV, 
and SPH MHD methods, implemented in the same code. This in¬ 
cludes problems such as the MRI and MHD jets which have been 
historically challenging. We show in all cases that the new nresh- 
free methods exhibit good convergence and stable behavior, and are 
able to capture all of the important behaviors, even at low resolu¬ 
tion. On some problem classes, we show they converge faster than 
state-of-the-art AMR codes using CT. 


2 NUMERICAL METHODOLOGY 

2.1 Review of the New Meshless Methods 

Paper I derives and describes the pure-hydrodynanric version of the 
numerical methods here in detail, including self-gravity and cos¬ 
mological integration. This is almost entirely identical in MHD; 


1 A public version of this code, including the full MHD implementa¬ 
tion used in this paper, is available at http : / / www .tapir . caltech. 
edu/~phopkins/Site/GIZMO. html Users are encouraged to mod¬ 
ify and extend the capabilities of this code; the development version of the 
code is available upon request from the author. 


therefore we will not repeat it. However we will very briefly review 
the new numerical methods. 

The equations we solve are the standard finite-volume 
Godunov-type equations: the fundamental difference between our 
meshless methods and a moving-mesh is simply that the definition 
of the volume partition (how the volume is divided among different 
nresh-generating points or “cells/particles”) is distinct. The further 
difference between this and a fixed-grid code is, of course, that the 
nresh-generating points/cells move, and that their arrangement can 
be irregular (as opposed to a Cartesian grid). 

In a frame moving with velocity Vframe, the homogeneous Euler 
equations in ideal MHD (and pure hydrodynamics, which forms 
the special case B = 0) can be written as a set of hyperbolic partial 
differential conservation equations of the form 

<9U 

+V-(F- Vframe® U) =S (1) 

where V ■ F refers to the inner product between the gradient oper¬ 
ator and tensor F, <g> is the outer product, U is the “state vector” of 
conserved (in the absence of sources) variables, the tensor F is the 
flux of conserved variables, and S is the vector of source terms 
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where p is mass density, e = u-\- |B| 2 /2p+ |v| z /2 is the total spe¬ 
cific energy ( u the internal energy), P T = P -\- |B| 2 /2 is the sum 
of thermal and magnetic pressures, and ip is a scalar field defined 
below. 

The meshless equations of motion are derived in Paper I in 
standard Galerkin fashion beginning from the integral form of the 
conservation laws, after multiplying Eq.|T]by a test function <j> 

0 = J dQ + J (F<j>)-n e nrfdfi (3) 

Here the domain f i is such that dQ = d^xdt, where v is the number 
of spatial dimensions, df /dt = df/dt + Vf rame (x, t) ■ V/ is the co¬ 
moving derivative of any function /, and ngo is the normal vector 
to the surface d D, and the test function is an arbitrary differentiable 
Lagrangian function. 

To transform this into a discrete set of equations, must chose 
how to partition the volume (for the “averaging/integration” step). 
If we choose a uniform Cartesian grid between uniformly spaced 
points, then we will recover the standard Godunov-type finite- 
volume grid-based equations of motion (like that in ATHENA and 
many popular AMR codes). If we choose a Voronoi tesselation be¬ 
tween moving mesh-generating points, we recover a moving-mesh 
method similar to AREPO. For the new methods in Paper I, we par¬ 
tition the volume according to a continuous weighting function /, 
such that the fraction of the differential volume d v x at the point x 
associated with the mesh-generating point at x = x,- is given by 


W(x-x h h(x)) 

J2jW{x-Xj,h(x)) 


Where W is any kernel/weight function and h is a “kernel length.” 
Note that this guarantees a “partition of unity” (the volume is per¬ 
fectly divided into cell-volumes Vi), and leads to a Voronoi-like par¬ 
tition, but with slightly smoothed boundaries between cells (which 
leads to some advantages, and some disadvantages, compared to 
moving meshes, where the boundaries are strict step functions). In 
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the limit where W goes to a delta function, the method becomes 
exactly a Voronoi-type moving-mesh method]^ 

This choice is combined with a second-order accurate moving- 
least squares matrix-based gradient operator, which has been uti¬ 
lized in many other methods (including grid-based codes; see 


Maron et al. 

2012, Tiwari & Kuhnertj 

2003; Liu et al. 2005| Luo 

et al. 20081 

^anson & Vila 2008a|"b| 

Mocz et al. 2014b). Eq. [3 


can then be expanded and analytically integrated to yield a second- 
order accurate set of discrete evolution equations: 


j t {VU), + 'E i V i j-A l j = (VS), (5) 


This is identical to the standard Godunov-type finite-volume equa¬ 
tions. The term ViU/ is simply the cell-volume integrated value of 
the conserved quantity to be carried with cell/particle i (e.g. the to¬ 
tal mass nti = Vi pi, momentum, or energy associated with the cell 
i); its time rate of change is given by the sum of the fluxes F, ; - 
into/out of an “effective face area” A;;, plus the volume-integrated 
source terms. The full mathematical derivation and expression for 
the A ij is given in Paper 1 (§ 2.1). 

We then use a standard MUSCL-Hancock type scheme for 
finite-volume Godunov methods to solve Eq.[5] This is commonly 
used in grid and moving-mesh codes l van^LeerJJ984 JToroJ1997 


Teyssier|2002| Fromang et al.|2006||Mignone et al.|2007||Cunning-| 

ham et ah|2 009[ Springel 2 010) ; it involves a slope-limited linear 

reconstruction of face-centered quantities from each mesh generat¬ 
ing point (cell “location”), with a first-order drift/predict step for 
evolution over half a timestep, and then the application of a Rie- 
mann solver to estimate the time-averaged inter-cell fluxes for the 
timestep. See Paper I(§ 2 and Appendices A & B) for details. The 
points then move with the center-of-mass gas velocity. 

In Paper I, we derive two variants of this method and im¬ 
plement them in GIZMO. First, the meshless finite-volume (MFV) 
method. This solves the Riemann problem between cells assum¬ 
ing the effective faces move with the mean cell velocity; this is 
analogous to a moving-mesh code, and includes mass fluxes be¬ 
tween cells. Second, the meshless finite-mass (MFM) method. This 
solves the Riemann problem assuming the face deforms in a fully 
Fagrangian fashion; in this case there are no mass fluxes. The two 
are formally identical up to a difference in the non-linear (second- 
order) error terms in the fluxes, provided the cells move with the 
gas velocity. In practice, each has some advantages and disadvan¬ 
tages, discussed below. 


2.2 Code Modifications for MHD 


Everything described above is identical in hydodynamics and 
MHD; and all details of the code (except those specifically de¬ 
scribed below) are unchanged from Paper I. 


2 Here and throughout this paper, we will define the kernel size A, at the 
location of cell/particle i as the effective cell side-length, based on the cell 
volume Vi. In 1D/2D/3D, this is h, = V), h, = (Vi/n) l l 2 , Ilf = (3 V//4tt) 1 / 3 , 
respectively. The volume V,- is calculated directly from the neighbor posi¬ 
tions defining the volume partition (see Paper I). This exactly reproduces the 
grid spacing if the particles are arranged in a Cartesian grid. Note that, in 
principle, VP(x, h ) can be non-zero at |x| > h. In MFM/MFV methods this 
has nothing to do with the effective cell/particle volume/size (conserved 
quantities are not “smoothed” over the kernel, but only averaged inside the 
a single cell of volume V) just like in a grid code), but instead reflects the 
size of the stencil (number of neighbor cells) between which fluxes are com¬ 
puted. As in grid codes, increasing the stencil size can increase diffusion, 
but does not directly alter the resolution scale. 


As usual for finite-volume Godunov schemes, we explic¬ 
itly evolve the conservative variables (VB), (integrated magnetic 
field over the volume partition corresponding to a mesh-generating 
point) and ( miji)i = f ptpdVj', primitive variables and gradients are 
then constructed from these (e.g. B, = (V B ), /V, as in Paper I). 

2.2.1 The Riemann Solver 


As in Paper I, we solve the ID, un-split Riemann problem in the 
rest-frame of the effective face between the cells. However we re¬ 
quire a Riemann solver that allows B^O. Since in the hydro case 
we use the HLLC solver, here we adopt the widely-used HLLD 
solver (Miyoshi & Kusano 2005}. This is accurate at the order 
required, and extremely well-tested (see e.g. Miyoshi & Kusano| 
|2005[ [Stone et al.|2008| and modem versions of e.g. RAMSES and 
ENZO). 

The frame motion is calculated for both the MFM and MFV 
methods as in Paper I. We emphasize that for our MFM method, it is 
required that we use a solver which explicitly includes the contact 
wave (i.e. the contact discontinuity, on either side of which mass 
is conserved). This is because in MFM we always solve the Rie¬ 
mann problem with frame moving exactly with the contact wave; 
attempting to simply find a frame where the mass flux vanishes in 
a simpler (such as the HFF or Rusanov) approximation leads to 
incorrect, and in many cases unphysical solutions. 

2.2.2 Signal Velocities and Time-Stepping 


As in Paper I, we limit timesteps with a standard local Courant- 
Fridrisch-Fevy (CFL) timestep criterion; in MHD we must replace 
the sound speed with the fast magnetosonic speed: 


AtcFL.i = 2 Ccfl 



max 
; • 
sig, 


MAX 


[v f , ij + Vfji — MIN ^0, (— 


(6) 

v,Hx,- Xj ) \-| 

| Xj Xj | )\ 

(7) 


(vt,u ) 2 % 


cl i+V%i+ \J{c 2 s i+V 2 A i) 2 -4cl ( V 2 A i (B; ■ Xij) 2 


(8) 


Here x,- ( = x, — X/ is the separation between two points, x is the unit 
vector x/|x|, c s ,i is the sound speed, va is the Alfven speed, hi is 
the effective cell size defined above, MAX, refers to the maximum 
over all interacting neighbors j of i, and v s i g is the signal velocity 
( |Whitehurst| 1995 [ |Monaghan| 1 997) . 

2.2.3 Divergence-Cleaning 


Ideally, this would complete the method; but the above method can¬ 
not ensure the divergence constraint. To do this, we must add the 
following source terms: 

S —Spowell T" SDedner (9) 


( ° ^ 


( 0 \ 

B 


0 

v B 

- 

B-(VV') 

y 



V o ) 


\ (V ■ B) pc\ + P't/’/r / 


The first term (Sp OW eii) represent the jPowell et ah] 41999) or “8- 
wave” cleaning, and subtracts the numerically unstable terms from 
non-zero V ■ B. This is necessary to ensure numerical stability and 
Galilean invariance - most problems will crash or converge to in¬ 
correct solutions without this. The second term (SDedner) follows 
the method of[Dedner et al.||2002), who introduce a conservative 
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scalar field ip which transports divergence away from the source 
and damps it. This is necessary to keep V ■ B low, minimizing the 
resulting errors. 

Following |Dedner et al.| ( |2002| and [Gaburov & Nitadori| 
< | 201 1 | >, it is straightforward to show that this leads to the follow¬ 
ing form for the discrete terms in our equation of motion (Eq. m 


( 


(VS), —(VV-B)r 


0 

B, 

v, B, 
v, 


\ 


V o 

/ 0 

0 

B,-(w^)r 

(Wi/,)* 

V (VV ■ B)* pi c^ j + (mip)i/Ti J 


(VV.B)*e-^4|A, ; 

j 


( 11 ) 

( 12 ) 

(13) 


where the B' and ip terms are defined below. 

We have some freedom to choose c*., and r,. Following|Tricco| 


& Price 


( 2012 ), we take Chj = o\j 2 v^f / 2 , where v^ g A f is the 


maximum signal velocity as described above ( 07 , is simply a con¬ 
venience parameter which re-normalizes the characteristic speed). 
The value v^f/2 is close to the fast magnetosonic wave speed, 
but also accounts for super-sonic cell approach velocities, which is 
critical for good behavior in highly supersonic compressions. We 
have experimented with variations in the dimensionless parameter 
(j/,, and find the best results for ot, = 1 (values 07 , <C 1 produce 
ineffective cleaning, values 1 lead to numerical instability). We 
take Ti = hi/(a p c Tt i) (where /?,- is the effective cell length defined 
above)|^For c T ,i (the damping speed), we have considered several 
choices, all of which give very similar results. These are detailed 
in Appendix [D] our default choice is Eq. |D 6 | which is closely re¬ 
lated to the local fastest possible signal velocity. We have also ex¬ 
perimented extensively with a p , and find a best compromise be¬ 
tween stability and diffusivity for values a p ~ 0.05 — 0.3; we adopt 
a p = 0.1 as our default in all problems here. 


Note that, unlike the original [Dedner et al.| ( |2002) formula¬ 
tion (and implementations in codes like AREPO and PLUTO), this 
means c*., and t, are spatially variable. This allows us to main¬ 
tain hierarchical timestepping (while forcing a constant a, imposes 
a global maximum timestep, a severe CPU cost penalty), and has 


many other advantages (see Appendix[B]l. But it is not a priori ob¬ 
vious that this will maintain stability. However, in Appendix |B] we 
discuss this in detail and derive a rigorous stability criterion, which 
should be satisfied by our choices above, designed so that c* and r 
are locally smooth (on the kernel scale). We confirm this stability 
in our numerical tests. 

We caution that the discrete source terms, particularly the 
Powell terms which subtract i-centered quantities, are not mani¬ 
festly antisymmetric between cell pairs ij. This means that mo¬ 
mentum and energy conservation in MHD are only accurate up to 
integration accuracy, times a term proportional to V ■ B (unlike in 
hydrodynamics, where conservation can be ensured at machine ac- 
curacy)f]Controlling V ■ B is critical to minimize these errors. 

These source terms also modify the Riemann problem. When 
we perform the reconstruction to obtain the left and right states 
Ul (J- side) and U„ (i-side), we can define a convenient coordinate 
system where x = A ,7 (i.e. the jr-axis is normal to the effective 
face between cells i and j). In this coordinate system, the normal- 
component of the 6 -fields B ' x , will in general not be equal. But equal 
values (i.e. non-zero V ■ B in the ID problem) are required for a 
physical solution. Without divergence-cleaning (Powell-only), this 
is handled by simply replacing B' x , and B' x R with the mean value 
B x jj = (B xL + B xR )/ 2.(Dedner et aL (2002) showed that with the 
source terms of Eq. [9] the infinitely-sharp discontinuity leads to a 
physical solution B x L —7 B' x R —7 6 ,., 7 , ipi. —7 ipR —7 1 / 1 , 7 , in infinites¬ 
imally small time: 




Ch, i 


2 

Vf,L 


B x ,i j — X (B X ,L + B x ,r) + —- (tpL — IpR) 


Ch, i 


= j (V’i + IpR) + 

= MAX [v f ,L, v f ,j?] 
1 2 2 

^,L + V A,L + 


2 Ch, ij 

(b x . l -b' r ) 


\J( c Il + ' a,l ) 2 - 4 c IlB x 2 l / p L 


(14) 

(15) 

(16) 
(17) 


here c*,y is the fastest wave speed in the local ID problem, which 
can be computed only from the i and j values (it does not necessar¬ 
ily correspond to c/ v in Eq. |9]above)|^]This is separable from the 
full Riemann solution. So, in the Riemann problem, we first update 
B x and ip according to the above, then compute the full Riemann 
solution using the updated values (and usual B' y LR and B' z LR ). The 
flux of B' x is then F tJ IP = v/, acc B' x l] , where Vj.face is the normal 
face velocity in the boosted frame (in which we solve the Riemann 
problem). Because ip is advected with the fluid, we follow |Gaburov| 
and simply take the ip flux to be F;y ^ = F t j, p ip R 


& Nitadori 


(2011 


for Fij, p > 0, and Fy,i/> = Fij, P ipR for Fj.p < 0, where Fy,p is the 

mass flux (so this vanishes for our MFM method)^] 

The HLLD solver requires an initial guess for the left and 


3 Eq .[9]is the continuum equation; we stress that we are not free to choose 
how we discretize it. To actually ensure numerical stability, the form of the 
V • B terms must exactly match those terms from our Riemann solver so¬ 
lution which are unstable (e.g. the tensile terms). Likewise, the divergence 
cleaning must act specifically to reduce V ■ B defined in the same manner, 
or it does not serve any useful purpose. It might be tempting, for example, 
to use the value of V • B calculated from our second-order accurate matrix 
gradient estimator for the Powell terms, or to construct a pair-wise sym¬ 
metrized version of Eq. |1 l|which manifestly maintains momentum conser¬ 
vation. However, these will not actually eliminate the unstable terms, and 
we have confirmed that they lead to catastrophic errors in our tests. 

4 In timestepping, we update (mi/)),- for the 77 term with the implicit solu¬ 
tion (mf))i oc exp (— Ati/ri); this allows us to take larger timesteps without 
numerical instability. 


5 Some of the Dedner terms can be made anti-symmetric without destroy¬ 
ing the numerical stability of the scheme; for example the fj-fiux correction 
described below, and the (VV • B) term (by using a single wavespeed 
c h.i = c h f° r the whole problem). However in all our tests the conserva¬ 
tion errors from those terms are always sub-dominant to the error from the 
(inescapable) Powell (V V • B)* v, term (and all these errors vanish when 
V • B —> 0). So we find the best overall conservation properties result from 
using the most accurate possible cleaning scheme, rather than a partially- 
conservative (but less accurate) cleaning. 

6 We have also explored an alternative, two-wavespeed formulation of the 
B and *0 terms, discussed in Appendix]?] For all tests here, the difference is 
small. 

7 We have experimented with using fjij for the ^-flux. However, this yields 
no improvement on any test problem, and the (B' x L — B' x R ) term from 
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Meshless MHD 5 


right wavespeeds, to compute a solution. If we define v' x as the 
normal-component of the reconstructed velocities, then we use 
Sl = MINfi^ z,. v' xR \ — MAX[vf L , v|- J K ], where v/ is the updated 
fast magnetosonic wavespeed using the updated normal 


(vjV) 2 = 


c s,L + V A,L + 


\f( C 


V A,lY - Ac lL B 'x 2 ij/pL 


(18) 


and Sr = MAX[< t , v' x , R ] + MAX[vjf t , v^][| 

Finally, as in Paper I, we solve the Riemann problem in the 
boosted frame Vframe corresponding to the mean motion of the 
quadrature point between mesh-generating points. We must there¬ 
fore boost back to the simulation frame. The de-boosted fluxes for 
cells follow Paper I, for the hydro terms, but with the additional 
terms for B (for the fluxes to cell-/): 


(F,y(B) • Ay) . ->fij (B) • Ay - B' xJj |Ay| Vframe (19) 

(Fy(e) • Ay). —>Fy(e) - Ay — B' xi j I Ay I Vframe • B; (20) 

The second equation just accounts for the energy flux associated 
with the corrected B-flux. 


2.3 The SPH MHD Implementation in GIZMO 


As described in Paper I, GIZMO is a multi-method code: users can 
run with the MFM or MFV hydrodynamic methods, or SPH, if 
desired. We therefore update our SPH implementation to include 
MHD. The exact SPH equations are given in Appendix [A] 

Briefly, the non-magnetic implementation of SPH follows the 
“modern” P-SPH method developed in Hopkins] ( ]2013| > and ex¬ 
tended in Paper I. This includes state-of-art re-formulations of the 
SPH hydrodynamics equations to eliminate the known “surface ten¬ 
sion” errors i jSaitoh & Makino|2013|[Hopkins|2013| , Lagrangian- 
derived terms to account for variable smoothing lengths ( [Springel] 
|& Hern quist 2002), addition of artificial diffusion terms for thermal 
energy ( [Price | | 2008] | Wads ley et al.| |2008), higher-order switches 
for artificial diffusion terms to minimize unnecessary dissipation 
l |Cullen &~D ehnen 2010), the use of higher-order kernel func¬ 
tions to allow larger SPH neighbor numbers and reduce the zeroth- 
order SPH errors fDehnen & Al y|2012 }, switches to prevent dis¬ 
parate time-stepping between neighbor particles ( | Saitoh & M akino 
2009), introduction of more accurate matrix-based gradient estima- 
tors ^Garcia-Senz et al.[2012) , and conservative, more accurate cou¬ 
pling of SPH to gravity ( |Price & Monaghan|2007||Bames|2012) . 

The SPH MHD implementation combines these improve¬ 
ments to SPH with the MHD algorithms from the series of pa¬ 
pers by |Tricco & Price| ( |2012| |2013| >; |Tricco| j2015). These in¬ 
troduce artificial diffusion for magnetic fields (artificial resistiv¬ 
ity), with a similar “switch” to reduce unnecessary dissipation, and 


can introduce numerical instability under some circumstances. Therefore 
we use the simpler ?/;-fl li x . 

8 The HLLD solver can fail in some very rare circumstances if “bad” 
guesses are used. We therefore check whether our first estimate produces a 
solution where the pressure is everywhere positive. If this fails, then we in¬ 
stead compute the Roe-averaged velocity vr oc and fast magnetosonic speed 
c Roe , and try Sl = MIN[v' L - v f j L , v Roe - c Roe ] and Sr = MAX[v' R + 
vj 1 R , v Roe + fRoel- We then check again; if this fails (which does not oc¬ 
cur in any test problem here) we test a Lax-Friedrich estimate (Sl = — Sr, 
with Sr from our first guess). If this somehow fails still, we go back and 
re-compute the interface using a piecewise-constant (first-order) approxi¬ 
mation, then check the series of wavespeeds again. If this fails, the code 
exits with an error (this only occurs when unphysical values are input into 
the solver). 



Figure 1. Linear magnetosonic wave test problem (see § |3.1) . Here a trav¬ 
eling, one-dimensional fast magnetosonic wave is propagated one wave¬ 
length; we then define the LI norm as the mean absolute error relative to 
the known analytic solution in density (top) or magnetic field B x ( bottom ; 
error here is equivalent to the numerical V • B 7 ^ 0 errors). We compare 
our new, meshless Lagrangian finite-volume Godunov methods (“meshless 
finite-mass” or MFM, and “meshless finite-volume” or MFV) from Paper I 
(see § |2.1) , to the best current implementation of SPH MHD (see § |2.3f , and 
to state of the art grid codes, here ATHENA run as a third-order PPM code 
using constrained transport (CT) to ensure V ■ B = 0 to machine precision. 
Dotted line shows second-order convergence (LI oc N~ 2 ); MFM/MFV and 
grid/PPM methods converge at this rate, as expected. Convergence is also 
good (LI oc N~ 23 ) in MFM/MFV for the divergence errors (( B x ) = 1, so 
these are fractionally very small). SPH shows some (slower) convergence 
until its known zeroth-order errors dominate; then errors flatten with resolu¬ 
tion. This is reduced in SPH by increasing the kernel size. “SPH-lo” (stan¬ 
dard A/ngb) uses the equivalent of Vngb = 32 in 3D (our default choice in 
all MFM/MFV runs). “SPH-hi” uses a 3D-equivalent /Vngb = 120. 

re-discretize the MHD equations from the particle Lagrangian in¬ 
cluded the |Dedner et al.|(2002) and |Powell et al.| < fl999) terms, so 
that the divergence cleaning actually acts on the tensile-unstable 
terms and these terms are properly subtracted (unlike most previ¬ 
ous SPH-MHD implementations). We make some further improve¬ 
ments: following |Price et aT]j2012) ; |Bate et al-H2014| and directly 
evolving the conserved quantities (VB),■ and ( mip)i , and slightly 
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6 Hopkins & Raives 


modifying the artificial resistivity terms to allow the method to cap¬ 
ture cosmological field growth and MHD fluid-mixing instabilities. 

3 TEST PROBLEMS 

We now consider a series of test problems. To make comparison 
as fair as possible, we will consider MFM/MFV/SPH implementa¬ 
tions in the same code (GIZMO). Unless otherwise stated, we com¬ 
pare to fixed-grid results from ATHENA (Stone et al. [2008fr ; this is 
representative of the state-of-the-art in finite-volume, non-moving 
mesh codes. Because we are interested in methods which can be 
applied to complicated, multi-physics systems, we use the same 
numerical implementation and identical values of purely numerical 
parameters (e.g. a p and Ccfl), within each method, for all prob¬ 
lems. It is of course possible to improve performance in any simple 
test problem by customizing/tweaking the method, but this usually 
entails a (sometimes serious) loss of accuracy on other problems. 
The implementations used here are therefore our attempts at a “best 
compromise” across all problems. 

3.1 Linear Magnetosonic Waves: Testing Convergence 

We begin by considering a simple linear one-dimensional magne¬ 
tosonic wave []The problem is trivial, but since virtually all numer¬ 
ical schemes become first-order at discontinuities, smooth linear 
problems with known analytic solutions are necessary to measure 
formal convergence. Following Stone et al. ( |2008| >, we initialize a 
unit-length (in the x-direction), periodic domain, with polytropic 
7 = 5/3 gas, density p = 1, pressure P = 1 /y, and magnetic field 
B/a/47t = (1. Vi, 1/2). We add to this a traveling fast magne¬ 
tosonic wavf^jwith amplitude Sp/p = 10 -6 , and allow the wave to 
propagate one wavelength; we then define the LI error norm for the 
density (or any other variable) Ll(p) = /V~* JA |p(jr,-, t ) — p(jq, t = 
0)|. We consider N = 16, 32, 64, 128, 256, 512, 1024, 2048. 

All methods we consider are able to evolve the wave. Fig. [T] 
plots the LI norm for density (the velocity variables look similar), 
and B x . In p, we see that both our MFM and MFV methods con¬ 
verge - as expected - with second-order accuracy. We compare 
to a state-of-the-art grid code, here ATHENA, run in the most ac¬ 
curate possible mode: PPM (formally a third-order reconstruction 
method), with the CTU integrator, and CT used to ensure V ■ B = 0. 
Despite the higher order of ATHENA and the fact that it uses CT to 
ensure V ■ B errors remain at the machine-error level, the errors 
are nearly identical to our MFM/MFV results. The LI norm for 
B x directly measures the divergence errors; since we do not use 
CT, these are non-zero. However they are (1) very small, and (2) 
converge away appropriately - in fact, we see super-convergence 

9 See http://www.astro.princeton.edu/~jstone/ 

Athena/tests/linear-waves/linear-waves.html 

For the grid and MFM/MFV methods, the perturbation is initialized by 
keeping the particles equidistant and modifying the conserved quantities. 
In SPH however, significantly better performance is obtained by initializ¬ 
ing the density perturbation by using exactly equal particle masses with 
perturbed locations (this allows apparent convergence in SPH to extend fur¬ 
ther). We therefore show this case for SPH (as it corresponds to the more 
likely case in real problems). We note that this initial condition also im¬ 
proves our MFM/MFV performance, but by a smaller amount. The SPH 
performance is also substantially improved if we “turn off” artificial vis¬ 
cosity (or set its minimum to zero), but this is numerically unstable and 
produces catastrophic errors in most of our tests (see e.g. ]Hopkins|2013| 
[MT5||Hu~eTal.|2014|[Rosswog|2014f. 


(LI oc N~ 23 ) for our MFM/MFV methods. As in Paper I, the con¬ 
vergence rate in all variables is independent of the kernel neighbor 
number in MFM/MFV: the choice only controls the normalization 
of the errors (larger neighbor numbers reduce noise, but increase 
diffusion). Based on our experiments in Paper I, we find roughly 
optimal results using a 3D-equivalent neighbor number /Vngb = 32 
(/Vngb ~ 4 in ID). 

For SPH, we see slower convergence in p, but it is reason¬ 
able at low resolution; but at high-resolution, the SPH zeroth-order 
errors (here, the “kernel bias” error in density estimation, and the 
systematic gradient error which results from imperfect particle or¬ 
der in the kernel) begin to dominate, and the errors flatten with 
resolution. This is more severe if we use a lower neighbor num¬ 
ber; going to higher-order kernels and higher neighbor numbers 
suppresses the errors, although they still eventually appear. For 
true convergence, /Vngb must increase with N, as is well-known 
(see |Zhu et aT||2014[ and references therein). Here we compare 
3D-equivalent /V N gb = 32 (our default MFM/MFV choice used in 
all runs in this paper), henceforth referred to as “SPH-lo,” to 3D 
/Vngb = 120 (henceforth “SPH-hi”)f3 

3.2 Brio-Wu Shocktube: Capturing MHD Discontinuities & 
Controlling Noise 

Next we consider the |Brio & Wu| (1988| > shocktube; this tests 
whether the code can accurately represent uniquely MHD shocks, 
rarefactions, and contact discontinuities. We initialize a 2D periodic 
box (size 0<x<4, 0<y< 0.25, with 896 x 56 cells/particles) 
with left-state (p, v x , v y , v z , B x , B y , B z , P) = (1,0,0,0,0.75,1,0,1), 
and right-state (0.125,0,0,0,0.75, —1,0,0.1), with 7 = 2. Figs. [2]- 
[3] compare the results at time t = 0.2. Note that our box is inten¬ 
tionally large and extends well beyond the “active” domain; the 
dynamically active region is only ~ 200 elements across □] 

At this resolution, the CT-based grid code is well-converged, 
except for some post-shock “ringing” most visible in v x . In MFM & 
MFV methods, the agreement with the exact solution is good at this 
resolution. At 4x higher-resolution, we find that the MFM/MFV re¬ 
sults are nearly indistinguishable from the exact solution. At lower 

11 Increasing the neighbor number /Vngb with particle number N must be 
done carefully in SPH, as discussed in e.g. |Price|(2012); |Dehnen & Aly| 
(20^ ; [ZhueTari(20T4t , since one must avoid the pairing instability and 
also account for the fact that simply increasing neighbor number in SPH 
changes the resolution length. For SPH-lo we use the popular [Schoenberg] 
(1946) cubic spline kernel, so /Vngb — 32 in 3D corresponds to an “effec¬ 
tive” resolution scale of h « 1.1 in units of the mean inter-particle spacing 
(followingjDehnen & Aly|2012| h = 2a, where a is the standard deviation 
of the kernel). For SPH-hi we use the [Schoenberg (1946]) quintic spline, 
so A/nqb — 120 in 3D corresponds to h « 1.4. Therefore we caution that 
the “effective resolution” of SPH-hi is slightly larger (« 20%) than SPH-lo 
at fixed N\ however the difference is much smaller than might naively be 
expected based on /Vngb alone. 

12 We have run this test with a number of different element configurations 
in our initial conditions, including: square, triangular, and hexagonal lat¬ 
tices, SPH and gravitational glasses, and random (Poisson) positions. We 
have also considered the case of equal particle masses (different particle 
spacing across the initial discontinuity) or unequal particle masses (same 
spacing). As expected from the derivation in |Gaburov & Nitadori| (20TT) , 
our MFM/MFV methods are only very weakly sensitive to these choices 
(slightly more noise appears in the “worst case” situation of random posi¬ 
tions). Perhaps more surprising, our SPH results are also only weakly sen¬ 
sitive to this. We show results for a relaxed SPH glass but they are qualita¬ 
tively identical to any of the other configurations. And we show explicitly in 
Paper I that the choice of equal/different particle masses makes a difference 
only in the magnitude of the errors just at the contact discontinuity. 


© 0000 RAS, MNRAS 000, 000-000 


































Meshless MHD 1 



- 1.0 - 


rN 

.. 





1.6 1.8 2.0 2.2 2.4 2.6 2.8 


2.0 - 


a 1.5 - 


1.0 — 


R 


1.6 1.8 2.0 2.2 2.4 2.6 2.f 



1.6 1.8 2.0 2.2 2.4 2.6 2.8 


1.6 1.8 2.0 2.2 2.4 2.6 2.8 


Figure 2. Brio-Wu shocktube (§ |3.2| , at time t = 0.2; we compare the exact solution to that computed at finite resolution (plotted region contains ~ 200 
elements across the x-direction) with different methods. High-order grid methods have converged well at this resolution, except for post-shock ringing in v x . 
MFM/MFV methods also show good convergence; but at this resolution, MFM still shows some small “overshoot” in the jumps atx ps 2.1, 2.3 (more sensitive 
to our slope-limiter than the method itself), and both show some small (percent-level) errors in B x owing to the V • B errors; however the fractional magnitude 
of V • B is controlled well by our cleaning scheme (typical errors ~ 10 —4 at this resolution; still below 10 -2 at jumps). Discontinuities are well-captured 
across ~ 2 cells/particles in the linear direction. SPH (with high Mmgb) captures all the key features, but at this resolution shows larger noise, ~ 20%-level 
overshoots in v x , v y ,B y , P at rarefactions ( x ps 1.9, 2.6), some suppression of internal energy around x ~ 2.3 (owing to more smeared-out dissipation from 
divergence-cleaning), and significantly larger V • B errors (reaching ~ 10%); these do converge away but more slowly. 




Figure 3. Brio-Wu shocktube, as Fig. [2] for additional methods. If we consider MFM/MFV using only the |Powell et al.|fl999) source terms to stabilize MHD, 
but no|Dedner et al. (2002} divergence-control (as has often been done in the literature), we obtain incorrect shock jumps; most noticeably in u & B x . This error 
does not converge away. This is despite the fact that the formal V • B errors are still small; the key is the terms in the |Dedner et al.|j2002) scheme that enter 
the Riemann problem and act specifically at discontinuities. We also compare SPH run with the same (lower) neighbor number as our MFM/MFV methods; 
here the noise is larger (as expected). 


resolution, there is some “overshoot” at the density discontinuity, 
with MFM - this is discussed in detail in Paper I; it is mostly sensi¬ 
tive to the choice of slope-limiter (not the basic numerical method). 
The effects are much smaller in MFV, owing to mass fluxes allow¬ 
ing more sharply-captured density discontinuity. This also causes 
a small pressure “blip” in MFM at the contact discontinuity, but 


this converges away. Shock jumps and discontinuities are captured 
across ~ 2 cells/particles in each direction; comparable to high- 
order grid methods. 

As an indicator of the V • B errors, we plot the dimensionless 
magnitude of hi |V B|//|B|,; our divergence cleaning keeps this 
generally at low values (<^ 10 -4 ), except at the magnetic shocks 
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8 Hopkins & Raives 



Figure 4. The Toth super-sonic shocktube (§ |3.3f The resolution is similar 
to the Brio-Wu test. MFM and grid/CT methods converge most rapidly to 
the exact solution (dotted), followed by MFV, which exhibits some residual 
noise in p around x ~ 2.1 — 2.4 at this resolution. MFM/MFV both control 

V ■ B errors well; the large B y discontinuity introduces a small (percents- 
level) offset in B x which converges away oc N~ 1 (nearly ideal), this is es¬ 
pecially challenging for divergence-cleaning methods. SPF1 shows larger 

V • B errors, and noise; with some systematic offset in the By jump around 
x ~ 2.2 (as divergence cleaning is spread over several smoothing lengths), 
which converges slowly. The noise is reduced with larger neighbor number, 
but not eliminated. In all methods, the Powell scheme alone leads to sys¬ 
tematically incorrect shock jumps in />', and [>: as in Fig. [2j these do not 
converge. 

(the large discontinuities in B y )\ but even there, the maximum value 
is still < 1(P 2 . The good divergence-cleaning is also manifest in 
B x ; at the same shocks, there are some jumps in B x generated, but 
this is returned to the correct value across ~ 2 particles, and the 
magnitude of the deviations from the analytic B x = 0.75 is typi¬ 
cally at the sub-percent level, at this resolution. 

With no divergence corrections whatsoever, the problem 


crashes. However, Fig. [3] shows that if we run with only the Pow¬ 
ell terms in Eq. [9] (i.e. do not include the |Dedner et al.| \2002) 
divergence-cleaning and damping terms), this particular problem 
is not badly corrupted in most respects. As expected the divergence 
and deviation in B x are larger. More seriously, though, a systemat¬ 
ically incorrect jump in u appears, which does not converge away 
without divergence cleaning, even at l(k higher resolution in the 
.r-direction. This problem occurs in MFM, MFV, and SPH. 

SPH captures all the qualitative features. However, if we run 
with the same neighbor number (kernel size) as in our MFM/MFV 
methods, the noise is larger. This is shown in Fig. [3] If we instead 
use the equivalent of a 3D neighbor number of ~ 120 (as opposed to 
the ~ 32 we use for MFM/MFV), the noise is reduced, as expected. 
However, there is still much larger noise and post-shock ringing 
(compared to MFM/MFV), and significant overshoot in the veloci¬ 
ties and u at the rarefactionsJ^Divergence-cleaning works in SPH, 
but is much less effective, especially around the contact discontinu¬ 
ity, where the divergence errors reach ~ 10%. We stress though that 
these errors are resolution-dependent, and do converge away even¬ 
tually. Better results (at a given resolution) can also be obtained on 
shocktube problems in SPH by increasing the strength of the artifi¬ 
cial dissipation terms (viscosity/resistivity/conductivity); however 
this significantly degrades performance on other tests we consider. 
For an example of this test which demonstrates good agreement be¬ 
tween SPH and the exact solution (by combining higher resolution 


and stronger dissipation), see Tricco & Price ( 

2013). 

We have also compared the Ryu & Jones 

(1995) MHD shock 


tube, which exhibits all seven MHD waves simultaneously. The 
qualitative results and differences between methods are the same 
as in the Brio-Wu test, so we do not show it here. 

3.3 Toth Shocktube: The Critical Need for Divergence 
Cleaning Beyond Powell 


Next we consider the |Toth| t [2000| l shocktube; this tests 
super-sonic MHD shocks. It is particularly important be¬ 
cause Mignone & Tzeferacos ( j2010) > showed that a hyper¬ 
bolic divergence-cleaning scheme such as the |Dedner et al.| 
( |2002| > method (or CT) is necessary to get the correct shock 
jump conditions at any resolution, in grid-based methods. 
We initialize a 2D periodic box with the same initial ele¬ 
ment configuration as the Brio & Wu (1988) 14 with left-state 


(p,v x ,v y ,v z ,B x ,B y ,B z ,P) = (1, 10,0,0,5/V47r,5/v^vr,0,20), 
and right-state (1, — 10,0,0,5/\/47r, S/a/T^O, 1) and 7 = 5/3. 

Fig. 0 compares the results at time t = 0.08 (for clarity, we 
plot only a randomly-chosen subset of 500 cells/particles for each 
method). Our MFM method does extremely well at this resolution; 


13 In extensive experiments, we found that the magnitude of these errors 
in SPH, at fixed resolution, is determined by the artificial viscosity and 
resistivity schemes. If we simply assume a constant (large) artificial vis¬ 
cosity (i.e. disable our normal “shock detection” switch) we are able to 
eliminate most of the noise, and reduce the overshoot, in SPH in Figs. [2} 
[3] However, such a choice would severely degrade the performance of our 
SPH implementation on almost every other problem we consider. As an 
attempted compromise, in AppendixfA] we discuss modifications to the de- 
fault |Cullen & Dehnen|(20l0) viscosity scheme and |Tricco & Price]j2Q13) 
resistivity scheme, which we have used here. Using instead the default ver¬ 
sion of the viscosity scheme makes only small differences for large neighbor 
number (“SPH-hi”), but produces much larger (order-unity) noise levels at 
low-neighbor number (“SPH-lo”) on this test. In highly super-sonic tests, 
the difference is negligible. 

14 Again, we have considered a variety of initial element configurations. In 
thejTdthj j2000j shocktube we find these produce negligible differences. 
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the only difference between it and high-resolution grid runs is a 
small (~ 1% level) deviation in B x introduced by V ■ B errors at this 
resolution, and some noise at the shock in the small v y . The V ■ B 
errors are generally extremely small, ~ 1% at super-sonic shocks 
and ~ 10 -4 elsewhere. MFV is similar, although there is substan¬ 
tially more noise at the same resolution in p and u in the post-shock 
region; this is seen in pure hydro in Paper I and in |Gaburov & Nita~2| 
|dori| ( |20lT| >; the larger noise translates to slightly larger V • B errors 
at shocks, hence slower convergence in B x . Both MFM & MFV are 
indistinguishable from the exact solution at lO.r larger resolution in 
the j-direction. 

Flowever, with the Powell-only (no divergence-cleaning) 
mode, we see that the shock is in the wrong place! The shock jump 
is systematically wrong in p, P and B y , and this leads to it being in 
the wrong place over time. Again, this appears even at infinite res¬ 
olution. The noise, especially in B x , is also much larger, and V ■ B 
errors, as expected, are factor ~ 100 larger. 

In SPH, the noise is much larger, even with much larger 3D- 
equivalent Angb = 120. This is especially noticeable in v y , B x , and 
B y . In B y , the shock jump is also systematically over-estimated by 
~ 2 — 5%, owing to the much larger (~ 10 — 30%) V ■ B errors at the 
shock jump. As discussed below (§|5j, in SPH, divergence-cleaning 
cannot act over smaller scales than a few smoothing lengths, so 
the method has difficulty controlling the errors seen in the Powell- 
only case. The divergence errors are systematically larger by factors 
~ 10 — 100. With (Vngb = 32 as in our MFM/MFV methods, the 
noise is yet larger (order-unity fractional noise in B x , v y ). 


3.4 Advection of a Field Loop: Minimizing Numerical 
Diffusion 

The next test is a standard test of advection errors and numer¬ 
ical dissipation. We initialize a periodic 2D domain: inside a 
circle of R = 'Jxr J ry 1 < Ro = 0.3 about the origin, we set 
(p, B x , B y ) = (2, Boy/R. — Bqx/R) with Bo = 10~ 3 . Outside the 
circle ( p,B x ,B y ) = (1,0,0), and everywhere (P,v x ,v y ,v z ,B z ) = 
(1,2,1/2,0,0); this is an equilibrium configuration that should 
simply be advected. In Fig. [5] we plot images of the magnetic en¬ 
ergy density. In Fig. [6] we plot the total magnetic energy in the box 
as a function of time (which should remain constant at its initial 
value). For the sake of direct comparison between grid and mesh- 
free codes, in these plots we take the initial element configuration 
to be a Cartesian grid, with unequal-mass elements. 

Advection of any configuration not perfectly aligned with the 
grid is challenging in grid codes; here the loop is continuously 
diffused away, at a rate that increases rapidly at lower resolution. 
In Lagrangian methods, on the other hand, stable configurations 
with bulk advection should be advected perfectly. In Paper I, we 
demonstrate that our MFM & MFV methods can advect arbitrary 
pressure-equilibrium hydrodynamic configurations (including arbi¬ 
trary scalar quantities) to within machine accuracy. However, here 
the introduction of the divergence-cleaning source terms leads to 
some initial diffusion of B in MFM/MFV. This is enhanced by the 
fact that the particle masses change discontinuously at the “edge” 
of the loop. But still, we clearly see the benefit of a Lagrangian 
method: convergence is much faster than in fixed-grid codes (even 
using CT); the dissipation in our 256 2 simulation with MFM/MFV 
is approximately equivalent to that in ATHENA at 1024 2 . In the im¬ 
age, we see slightly more noise; this is the expected “grid noise” 
which is higher in meshless methods, but the diffusion is less (in 


Initial Conditions Grid 




0 0.25 0.5^JV7^^1 

MFV MFM 




SPH-hi Powell 



Figure 5. Field loop advection test (§ |3.4| . An equilibrium, 2D field loop 
is advected uniformly across the domain; we plot |B| 2 (in units of the ini¬ 
tial loop value, as labeled) at time t = 20, for tests with 256 2 resolution. 
Here we use in initially Cartesian particle lattice with unequal-mass parti¬ 
cles (the density within the loop is a factor ~ 2 higher than the background). 
This is the natural configuration for fixed-mesh codes but represents a near 
“worst case” for Lagrangian codes. In all cases, the initial conditions (top 
left) should be reproduced identically. In non-moving grid methods, even 
at arbitrarily high order, advection errors diffuse the loop, while numeri¬ 
cal resistivity reduces the central field strength. In MFM/MFV/SPH, the 
advection errors are eliminated, and numerical resistivity at the center is 
reduced; however divergence-cleaning and “grid noise” from different par¬ 
ticle masses around the loop edge produce some diffusion and noise (the 
peak value of /z|V • B|/|B| in any cell at any time remains <C 0.01, how¬ 
ever). The Powell-only scheme exhibits much more severe noise, because 
the V • B errors are transported but not damped; this leads to non-linear 
corruption of the solution. 

particular, the “hole” which app ears at the center owing to numeri¬ 
cal resistivity is minimized) rj 


15 Of course, in any of our Lagrangian methods, it is trivial to obtain per¬ 
fect evolution of the field loop test for arbitrarily long times, if we simply 
disable any dissipation terms. In MFM/MFV this amounts to invoking the 
“energy-entropy” switch described in Paper I (setting it always to “entropy,” 
i.e. using purely adiabatic fluxes), and in SPH it amounts to disabling the ar¬ 
tificial dissipation terms (see examples in e.g. Rosswog & Price 2007 Price 
|2012) . Then, because the system is in uniform motion there is zero advec¬ 
tion and the evolution is trivial. The same is true for moving-mesh codes. 
But these changes make the methods numerically unstable (and would pro¬ 
duce disastrous errors at any shocks or discontinuities). We therefore will 
not consider such modifications further. 
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Figure 6. Quantitative decay of the box-averaged magnetic energy in the 
field loop test (Fig. |5J, owing to numerical diffusion/resistivity. At infi¬ 
nite resolution, methods preserve the initial (|B| 2 ). Because they are La- 
grangian, MFM/MFV/SPH methods show much less dissipation than high- 
order grid methods at the same resolution (default runs here are 256 2 , but 
we compare 64 2 — 512 2 MFM/grid runs for reference). SPH shows some 
spurious initial growth of (|B|) as divergence-cleaning acts on zeroth-order 
kernel errors, but the subsequent decay rate is close to MFM/MFV. Other¬ 
wise on this test SPF1 is not as sensitive to (Vngb ■ Powell-only methods are 
unstable on this problem, and lead to artificial field amplification. 


If we use Powell-only divergence subtraction, the conserva¬ 
tion errors associated with non-zero V ■ B can actually lead to non¬ 
linear growth of B, such that the total magnetic energy increases! 
The growth in this case is nearly resolution-independent, since it is 
sourced around the sharp discontinuity in the field at the edge of the 
loop. By the end of the simulation, the total magnetic energy in the 
Powell-only run has increased ~ 50%. The noise and asymmetry 
in the image have grown severely. 

In SPH, there is an initial, brief but unphysical growth in 
|B| 2 ; this comes from the divergence-cleaning being less effective 
than in MFM/MFV (so it behaves like the Powell case); however 
once enough diffusion and particle re-arrangement has occurred, 
the divergence-cleaning operator can work effectively, and the en¬ 
ergy decays at approximately the same rate as our MFM/MFV cal¬ 
culations. 

As noted above, the noise and dissipation in the mesh-free 
methods is enhanced by the artificial jump in element masses at 
the edge of the loop (this is not a “natural” configuration for our 
mesh-free methods). We therefore consider in Fig. [7]the same test, 
with the same total number of elements in the box, with an initial 
triangular lattice configuration and equal-mass elements. The com¬ 
bination of reduced noise at the loop edge, and a factor ~ 1.25“ 
higher resolution within the loop (since it is higher-density), does 
reduce the noise and errors significantly. However, the qualitative 
behavior in every case is identical. 




Time 

Figure 7. As Figs. |5|6| but with a different initial particle configuration: 
here a triangular lattice with constant particle mass (with fixed total particle 
number, so the resolution is a factor of (\/2) 2 higher within the magnetized 
loop). In this case, the errors are reduced in all our Lagrangian methods. 
In particular, the noise in SPH, which is sensitive to both the local particle 
arrangement and differences in particle masses, is greatly reduced. Qualita¬ 
tively, the features in all cases are identical, however. 


3.5 Hawley-Stone Current Sheet: Numerical Stability 

This test follows |Hawley & Stone| < [1995| l. In a 2D periodic 
domain with —0.5 < x < 0.5, —0.5 < y < 0.5, we initialize 
(p, P, v x , v y , Vj, B x , B z ) = (1, P/2, Asin(27ry), 0, 0, 0, 0) and 7 = 
5/3, with B y /{AitY^ 2 = 1 for |jc| > 0.25 and B y /{ATt) 1 ^ 2 = —1 oth¬ 
erwise. This is not a good test of algorithm accuracy, since the non¬ 
linear solution depends sensitively on the numerical dissipation in 
different methods. However, it is a powerful test of code robust¬ 
ness. Qualitatively, the solution should exhibit rapid reconnection 
along the initial current sheet, which will launch nonlinear polar¬ 
ized Alfven waves, that generate magnetosonic waves, while mag¬ 
netic islands form, grow, and merge. For smaller /? and larger A, 
it becomes more difficult for algorithms to evolve without crash¬ 
ing or returning unphysical solutions (e.g. negative pressures in the 
Riemann problem). 

Fig.[ 8 ]shows the magnetic topology at time t = 5 in a run with 
/3 = 0.1, A = 0.1. For these parameters, MFM, MFV, SPH, and 
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Figure 8. Hawley-Stone current sheet (§ |3.5 . We plot magnetic field lines 
(arrows indicate local field strength and direction) at time t = 5 with 0 = 
0.1. A = 0.1, in MFM. MFV, SPH, and grid-methods produce very similar 
results for these parameters. Reconnection along the current sheets leads 
to magnetic “islands” which grow and merge. Our new mesh-free methods 
are able to stably evolve the current sheet indefinitely, with h |V ■ B|/|B| <S 
10~- for all cells at all times. 
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Figure 9. Stability limits of the current sheet problem in Fig. [8] Given an 
initial pressure P = (8/2, and velocity perturbation with amplitude A, we 
consider the minimum 0 and maximum A for which the problem can be 
evolved stably to time t = 10 (further to the top-right is more-stable). In 
the grid-based CT method of ATHENA, the total-energy formulation of the 
code, coupled with high-accuracy subtraction needed for accurate CT, and 
advection errors when the fluid moves over the grid, mean that the method 
will crash (negative pressures result) for 0 < 0.01 or A > 3. Combining the 
duel-energy formalism from Paper I, with a Lagrangian method that moves 
with the fluid, and using divergence-cleaning instead of CT, we are able to 
stably evolve the system until 0 reaches machine-error levels ~ 10 “ 16 , and 
arbitrarily large A > 10 5 (we have not considered larger A only because the 
simulations become too expensive, not because they crash). 



Figure 10. The Orszag-Tang vortex (§ |3.6| . We show images of density p 
(in code units, as labeled) at time t = 0.5, in runs with 256- elements (parti¬ 
cles/cells). All methods develop the major qualitative features, though there 
is some additional smoothing in SPH. Note that the contact discontinuities 
and shocks are captured sharply in MFM/MFV methods. 


grid methods (here, ATHENA) all look very similar^The real test 
arises when we vary 0 and A; in Fig. [9] we plot the maximum A 
and minimum 0 which we are able to use in each algorithm before 
the code crashes or returns an unphysical result. ATHENA crashes 
after some small early-time evolution for 0 < 0.01 or A > 3. The 
low-/? problem most likely owes to the fact that the method evolves 
total energy: when the magnetic energy dominates, this causes se¬ 
rious difficulty recovering the correct internal energy (since we 
must subtract two large numbers), eventually producing negative 
temperatures. Here we use a dual-energy formalism described in 
Paper I. which does not conserve total energy to machine error, 
but can handle essentially arbitrary ratios. So for our SPH, MFM, 
MFV implementations, /3 m j n is limited by essentially machine error 
(/3min ~ 1(T 14 — 10 -ls , depending on the formulation). For A, we 
find similar results; the increased stability owes both to the same 
dual-energy formalism above, but also to the Lagrangian method, 
which eliminates the advection errors that, in grid-based codes, be¬ 
come larger with the local fluid velocity. In non-moving grid codes, 
eventually, at any resolution, there is some bulk velocity which 
will wipe out the correct physical solution completely (necessitat¬ 
ing still-higher resolution); this is avoided in Lagrangian methods. 
We explore only values of A up to ~ 10 s in this test because the 
timestep becomes so small that it is impractical to evolve the sys¬ 
tem to late non-linear times, but we suspect the robustness of the 
algorithms should hold to similar machine error levels (i.e. allow¬ 
ing A ~ 10 16 ). 


16 For extensive description of this test problem in ATHENA, see http: 
//www. astro.Virginia. edu/VITA/ATHENA/cs .html 
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Figure 11. Comparison of the Orszag-Tang problem from Fig. |1()| We 
plot density p, magnetic field components B x , B y , and V • B, in horizon¬ 
tal slices at yo — 0.3125. With 256 2 cells, MFM, MFV, and grid methods 
have converged well to the exact solution (except for some small smooth¬ 
ing of the sharpest features). Here, the “exact” line is a 2048 2 ATHENA 
PPM CT result; at resolution > 512 2 , MFM/MFV/unboosted-ATHENA re¬ 
sults are indistinguishable. SPH performs well but shows further smoothing 
which converges more slowly. We also consider a grid simulation where 
the fluid is given an additional boost (a uniform v x = 10); the Lagrangian 
(MFM/MFV/SPH) methods are invariant to these boosts. But in grid codes 
the boost produces a smoothing of the B features at x ~ 0 — 0.2; these re¬ 
quire grid resolution > 1024 2 to converge away. Using only Powell diver¬ 
gence subtraction leads to a systematic error in p which is small, but does 
not converge away. 


3.6 Orszag-Tang Vortex: Shock-Capturing & Super-Sonic 
MHD Turbulence 

The Orszag-Tang vortex is a standard MHD test which captures a 
variety of MHD discontinuities, and develops super-sonic MHD 
turbulence, which is particularly challenging for many methods. 
In a periodic 2D domain of unit size, we take 7 = 5/3 and set 

(p, P, v,, v y , v z , B x , By, B z )=( 25/(36tt), 5/(12tt), - sin (2 ttv), 
sin(27rx), — sin(27ry)/(47r) 1 ^ 2 , sin(47^.x')/(47^) 1 ' ,2 , 0). 


Fig- [10] 
t = 0.5; Fig. 


compares images of the resulting density field at 
TT| quantitatively compares these by plotting the den¬ 



0.1 0.2 0.3 0.4 


Figure 12. Comparison of the Orszag-Tang problem as Fig. |10| but at time 
t = 1, when parts of the flow have broken up into chaotic turbulence. Unsur¬ 
prisingly the differences between methods have grown, but the results are 
still qualitatively similar. Note that the additional SPH smoothing in Fig. | 10| 
has here suppressed the formation of the central compact vortex; this feature 
is most sharply resolved in our MFM/MFV runs. 


sity and x, y components of B in a slice through y = 0.3125 at the 
same time. At this (early) time, the flow contains a complicated set 
of shocks, but has not yet broken up into turbulence. Figs 1 1 2| 13] 
consider the same at f = 1, when parts of the flow begin to evolve 
chaotically (amplifying the differences between methods). All runs 
here have 256 2 resolution. In all the methods here, all of the key 
qualitative features are captured at t = 0.5, and most at t = 1, in¬ 
cluding several complicated shocks, discontinuities, and sharp fea¬ 
tures. At t = 0.5, non-moving grid, MFM, and MFV solutions are 
essentially indistinguishable they have converged very well to the 
exact solution already, with only a slight smoothing of the sharpest 
features (as expected). Their “effective resolution” appears identi¬ 
cal. At t = 1, qualitative features are similar but quantitative dif¬ 
ferences appear; the MFM/MFV methods better resolve the central 
vortex/density peak and sharpest features in B (not surprising, since 
their Lagrangian nature automatically provides higher resolution in 
this region), but there are stronger oscillations in the low-density 
regions. However, if we boost the system by a constant velocity, at 
fixed resolution the stationary-grid solution degrades; for a v x = 10 
boost it is more comparable to a 64 2 run at both times. Obviously 
the Lagrangian methods avoid this source of error. Divergence er¬ 
rors are well-controlled at both times even in super-sonic shocks. 

SPH does nearly as well, although the features in B are more 
smoothed (similar to a 100 2 MFM/MFV/grid result) at t = 0.5, and 
this suppresses the appearance of the central density peak at t — 1. 
On this problem, there is not a dramatic difference, interestingly, 
between SPH with normal versus large neighbor numbers. As we 
saw in the shocktube tests, the divergence errors are larger by a 
factor of several in SPH compared to MFM/MFV. 

If we use a Powell-only cleaning scheme, the results are not 
too badly corrupted by V ■ B errors; however we do see some sys¬ 
tematic offsets, particularly in p around v ~ 0.6, and the position of 
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Figure 13. Quantitative comparison of the Orszag-Tang problem as Fig. ED 
but at time f = 1. Here we take the slice through the center (y = 1 /2), where 
differences are maximized. The density plot demonstrates the presence of 
the central vortex in particular; it is suppressed by smoothing in the grid 
code and disappears in the boosted grid, SPH, and Powell results at this 
resolution. 


the B jumps around x « 0. 1 at t = 0.5, and the central density peak 
and B x discontinuity at t = 1, that do not appear to converge away. 

Fig. |T4] demonstrates the convergence in our MFM method 
(MFV results are nearly identical here). As we increase the resolu¬ 
tion, we clearly see good convergence towards the exact solution in 
all features here. Quantitatively, the convergence in the Lj and Li 
norms of the plotted quantities is first-order, as expected due to the 
presence of shocks (the same is true in ATHENA). Compare this to 
older SPF1 MF1D implementations, which only saw convergence in 
some features, while others converged slowly or not at all (a well- 
known issue in SPFI; Stasyszyn et al.|2013) . 

We have also considered the pure driven-turbulent box prob¬ 
lems from Paper I (rms Mach numbers ~ 0.3 and ~ 8 ), as well as 
the ABC dynamo ( [Arnold et al.|198l) with small initial seed fields; 
we confirm that the turbulent power spectra and growth rate of mag¬ 
netic energy in the box agree well between MFM, MFV, and grid 
methods. This echoes our conclusions for the pure hydrodynamic 
case in Paper I. With Powell-only cleaning, however, the growth 


rate of the magnetic energy in the highly super-sonic case is artifi¬ 
cially high (growing faster than the flow-crossing time, indicating 
a clear numerical artifact). In SPH, the results depend on neighbor 
number, in a manner that we demonstrate in detail for the MRI test 
problem below. 

3.7 The Magnetic Rotor: Torsional MHD Waves 

Next we consider the MHD rotor, a standard problem from B 
|sara & Spicer| (1999 ) used to test strong torsional Alfven waves. 
A 2D domain of unit size and 7 = 7/5, is initialized with uniform 
P = 1, B = ( 5 /\/ 47 r, 0, 0). Inside a circle of size Ro = 0.1, p — 10 
and v = (—2y/Ro, 2x/Ro , 0); this is surrounded by a ring-shaped 
transition region from Ro < R < Ri = 0.115 with p = 1 +9 f(R), 
v = (-2 yf(R)/R, 2xf(R)/R, 0) with/(R) = (R, -R)/(Ri-Ro). 
At R > Ri, p = 1 and v = (0, 0, 0). 

As with the Orszag-Tang vortex, Figs. 1 15|T6] plot images of 
the magnetic energy density and field values in horizontal slices, at 
time t = 0.15. Again, all methods capture the key qualitative fea¬ 
tures. MFM/MFV and grid methods are very similar. Grid meth¬ 
ods converge slightly more rapidly on the sharp B field “spikes” at 
x = 0.3, 0.7 if the bulk velocity is nil, but if the fluid is boosted 
appreciably, the density spikes are noticeably more smoothed. In 
either case, MFM, MFV, and grid methods exhibit convergence at 
the same order. Divergence errors are again well-controlled, even 
around large discontinuities in B. 

Here, for SPH, using the same neighbor number as 
MFM/MFV (3D-equivalent 32 neighbors) leads to significant 
noise, and some systematic errors in B around x ~ 0.15 — 0.3, 0.7 — 
0.85. These are resolved if a large neighbor number is used. How¬ 
ever in both cases significant smoothing of the local extremum 
in the magnetic pressure (the brown patches at y ~ 0.5 ± 0.15 in 
Fig . 1 1 5[ > is evident, and V ■ B is less accurately suppressed. 

With only Powell-cleaning, significantly more noise is evident 
in the image; moreover, the systematic offset of the discontinuity 
positions in B x is clear, and this does not converge away. 

3.8 Magnetized Blastwave: Strong Shocks & 

Grid-Alignment Effects 

Next we consider a strong blastwave in a magnetic field; this is an¬ 
other standard problem, which tests the ability of the code to handle 
strong shocks and preserve symmetry. In the hydrodynamic version 
of this test, it is also very challenging for grid-based codes to avoid 
strong grid-alignment and preferential propagation of shocks along 
the grid (see Paper I); however these effects are reduced by the 
asymmetric forces in the MHD problem. 

We initialize a 2D periodic domain of unit size with 7 = 5/3, 
zero velocity, p = 1, B = (1 /x/2, 1 /x/2, 0), and pressure P = 10 
within a circle at the center of initial size R < Ro = 0.1 and P = 0.1 
outside. Figs. |17|l8| show images and slices through the solution at 
time t = 0.2, for 256 2 runs. Visually, the MFM/MFV, SPH (with 
high neighbor number), and grid (with zero boost to the fluid) so¬ 
lutions look similar. The blast expands correctly, with the cavity 
growing more rapidly in the direction along the field lines. There 
is some additional detail visible in the high-density regions in the 
Lagrangian solutions (the “dimpling” in the upper-right and lower- 
left); this is real and appears at slightly higher resolution in the grid 
calculation as well. The SPH solution is slightly more smoothed 
along the shock fronts. Here we also show the Powell-only and 
boosted-grid results; unlike the Orszag-Tang and rotor problems, 
here the differences are plainly visible by-eye. The Powell case de¬ 
velops incorrect features; the boosted-grid case loses symmetry. 

Quantitatively, we see these effects in Fig. [18] Note that 
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Figure 14. Resolution study of the Orszag-Tang vortex in our MFM method (MFV is essentially identical). We show both images as Fig. [lO] and slices as 
Fig. im for several resolutions. Most major features are present even at low resolution; convergence with higher resolution is clear for all features. The formal 
L\ and Lt convergence accuracy is close to ideal scaling oc IV -1 for this problem (given that it includes shocks and discontinuities). 
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Figure 15. Magnetic rotor (§ |3.7) . We show images as lie. 1 1f) | oI the mag¬ 
netic pressure |Bp/2 (in code units, as labeled), for runs with 256 2 resolu¬ 
tion. Most features appear identical; however note some difficulty in SPH 
capturing the pressure extrema at y ps 0.5 it 0.15, and additional noise in 
SPH (especially with low neighbor number). With Powell-only cleaning, 
similar errors and noise appear. 

MFM/MFV methods exhibit very small V ■ B values; SPF1 exhibits 
larger V ■ B, but still quite small. In general, different methods agree 
fairly well quantitatively. Flowever, for the Powell-only case, the 
value of B x , in particular, is seriously wrong in the upper-right por¬ 
tion of the solution. We emphasize that the error is worse if we re- 
simulate this with Powell-only cleaning but a resolution of 1024 2 . 
The failure of this method to guarantee the correct shock jumps cor¬ 
rupts the entire late-time solution. Note that, in SPFI (using the full 
|Dedner et al.| ( j2002j l cleaning), the magnitude of the V ■ B errors 
is not much smaller than our Powell-only run, but this incorrect 


behavior does not appear. This emphasizes that the magnitude of 
h |V-B|/|B| is not, by itself, the whole story; non-linear errors are 
damped by the |Dedner et al.| ( |2002[ > approach (at the cost of some 
additional local diffusion) which would otherwise build up coher¬ 
ently at later times. Even quite large values of V ■ B ~ 0.1 — 1 can 
be tolerated, if the algorithm includes a proper damping formula¬ 
tion (as our MFM/MFV/SPH implementations do). 

3.9 MHD Rayleigh-Taylor & Kelvin-Helmholtz Instabilities: 

Fluid Mixing in MHD 

Fluid mixing instabilities are astrophysically important and histor¬ 
ically challenging for SPH methods; the hydrodynamic forms of 
these are discussed at length in Paper I. In MHD. non-zero |B 
suppresses the growth of small-scale modes. If magnetic fields are 
too strong, no interesting modes grow. If fields are too weak, the 
case is essentially hydrodynamic. But there is an interesting, MHD- 
specific regime when the fields strengths are near-marginal; we 
consider this in a Rayleigh-Taylor (RT) problem. 

We take initial conditions from Paper I and |Abel| ( |201 l) l: 
in a 2D domain with 0 < x < 1/2 (periodic boundaries) and 
0 < y < 1 (reflecting boundaries), we take 7 = 1.4, B/t/ 47 t = 
(0.07, 0, 0), and p(y) = pi + [pi - pi)/(l +exp[-(y-0.5)/A]) 
with pi = 1, p 2 = 2, A = 0.025. Initial pressures are assigned to 
produce a gradient in hydrostatic equilibrium with a uniform grav¬ 
itational acceleration g = — 1/2 in the y direction (at the inter¬ 
face, P = P2/7 = 10/7 so c s = 1). An initial y-velocity perturba¬ 
tion v y = <5y v (1 +cos(87r(jf+ 1/4))) (1 +cos (57 t (y— 1/2))) with 
Sv y = 0.025 is applied in the range 0.3 < y < 0.7 (otherwise the 
velocities are zero)[] 

Fig. 1 19| shows the resulting density field at intermediate and 
late times, in a 128 x 256 simulation. In MFM, MFV, grid/AMR 


17 Following the method in Paper I, we use the routines generously pro¬ 
vided by R. O’Leary to construct the mesh-free initial conditions, then re¬ 
interpolate this onto the ATHENA grid. This is critical to ensure that the 
same “seed” grid noise is present in both methods, which in turn is neces¬ 
sary to see similar behavior in the late-time, non-linear phase of evolution. 
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Figure 16. Density and magnetic fields in a slice through yo = 1/2 in the 
MHD rotor problem in Fig. | 15 1 The qualitative differences between methods 
are the same as in Fig. |1 1 1 (the “grid+boost” run uses the same bulk v x = 
10). MFM/MFV/grid methods agree well. Note that boosted grid methods 
smooth the pressure spikes significantly, and Powell-only cleaning leads to 
systematically incorrect shock positions. 


with non-moving fluid, and SPH (with sufficiently large neigh¬ 
bor number), the linear growth of the field is essentially identical; 
this is consistent with our pure-hydro results in Paper I. Even the 
non-linear, late-time results agree reasonably well (although there 
is some symmetry-breaking in SPH which owes to less-accurate 
V ■ B-cleaning, even at large neighbor number). There is slightly 
more “grid noise” in MFV at this resolution (which leads to small 
differences in the late-time evolution). Just like in the pure hydro- 
dynamic case, in SPH, the results are very sensitive to neighbor 
number: if a smaller neighbor number is used (say, the same as we 
use in our MFM/MFV methods), then the initial seed mode is too 
weak - it is overwhelmed by the zeroth-order errors in the method, 
and no modes grow. Similarly, if we boost the fluid by a constant 
velocity in stationary-grid methods, advection errors break symme¬ 
try and dramatically suppress the growth of the mode at this reso¬ 
lution (much higher resolution is required to match the accuracy of 
the Lagrangian methods). All of this is consistent with our results 
from the hydrodynamic case in Paper I. 

If we consider the Powell-only case, the linear mode evolu¬ 
tion appears reasonable - the growth rates are only slightly sup¬ 
pressed. However, in the non-linear phase, the solution is totally 



Figure 17. Blastwave in a strongly magnetized medium (§ |3.8| . We plot 
density p (as labeled) in 256 2 simulations. The blastwave expands asym¬ 
metrically along field lines (initial B = (1/\/2, 1 / v2, 0)). as expected. All 
methods capture the key behaviors; Lagrangian methods capture slightly 
more detail in the high-density regions, but with enhanced grid noise. 
Stationary-grid methods converge more slowly when the fluid is boosted 
(here, v* = 10); the errors break the symmetry of the solution noticeably at 
this resolution, but they do eventually converge away. Using only Powell/8- 
wave methods (no V • B-damping) leads to physically incorrect shapes of 
the high-density features in the top-right and bottom-left comers; these do 
not converge away with resolution. 


corrupted! The non-linear errors accumulate, if only Powell-type 
schemes are used (once again, because in this method, divergence 
errors are only transported, not suppressed), until they overwhelm 
the real solution. Clearly, tests restricted to the linear regime are not 
sufficient to validate divergence-cleaning schemes. 

In SPH, we also find another difficulty unique to MHD. Nu¬ 
merical stability in SPH requires somewhat ad-hoc artificial dissi¬ 
pation terms for B, the “artificial resistivity.” As discussed in Ap¬ 
pendix El and |Rosswog & Price| ( |2007) ; |Tricco & Price| ( |2013| ), 
this carries ambiguities that are not present in hydrodynamics: the 
appropriate “signal velocity” for the resistivity could be the sound 
speed, Alfven velocity, or magnetosonic speed, and in some cases 
resistivity should be applied in rarefactions. The correct answer de¬ 
pends on the type of MHD discontinuity. In Godunov methods such 
as MFM/MFV/grid codes, the correct form of the dissipation is pro¬ 
vided by the appropriate solution to the Riemann problem. But in 
SPH, even with the state-of-the-art “switches” used here, this is 
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Figure 18. Slices through the MHD blastwave in Fig. m Most methods 
agree well; the convergence is good at this resolution except for SPH with¬ 
out a large neighbor number, or grids if the fluid is moving. The errors from 
the Powell-only (no-cleaning) method are dramatic: in B x , we see factor 
~ 2 — 3 systematically incorrect results. 

difficult to assign correctly. By default, we use the Alfven velocity 
in this switch. But |Tricco & Price] \2i) 1 show that this can lead 
to too-low a resistivity in super-sonic MHD turbulence, in turn pro¬ 
ducing shock breakup and serious noise (see their Fig. 7). They rec¬ 
ommend increasing the dissipation by using the fast magnetosonic 
speed Vfast instead. However, if we do that for this problem, it pro¬ 
duces excessive dissipation of the magnetic field around the contact 
discontinuity^] This makes the problem behave (incorrectly) like 
the pure-hydrodynamic case. 

In Fig. [20] we show full 2D maps of the divergence h |V ■ 
B|/|B|. in MFM/MFV, these are extremely well-controlled, with 
median values < 10 -4 and maxima < 10 -2 . In our default SPH 
implementation, they are larger by a factor of a few. In the “vf ast 
AR” SPH run, we see much larger values along the contact dis¬ 
continuity. These are not, however, caused by poor V ■ B-control; 
rather, excessive dissipation of B around the discontinuity leads to 


18 For the test in Fig. |l9| labeled “Vf ast AR,” we keep everything about our 
SPH MHD method identical (including the “switch” and maximum viscos¬ 
ity olb = 0.1), except replace the Alfven velocity in Eq. |A9| with the fast 
magnetosonic speed. 


a local sharp depression of |B| (the denominator), as opposed to 
| V - B|. In the Powell runs the divergence errors are (as expected) 
much larger. 

We have also compared the magnetized Kelvin-Helmholtz in- 
stabilit y, sho wn in Fig. [2T] The initial conditions follow jMcNallyj 


et al. 


: y, sh own in Fig. |zij 1 lie initial conditions tollow McJNaliy| 
(2012' (see also Paper I), a 2D setup with 256 2 resolu¬ 


tion elements following two streams with P = 5/2, 7 = 5/3, and 
(p, v x ) = (1,—0.5) and (2, 0.5), with a 1% amplitude initial seed 
mode and small interface region between the streams. We add a 
uniform B = (0.1,0, 0), about a factor « 2 below the critical value 
which suppresses the instability, and show results at f = 1.6 and 
t = 3.2 (where the KH growth time is « 0.7). Here, we obtain qual¬ 
itatively identical conclusions to our RT test. In the linear and even 
non-linear stages, MFM/MFV and non-boosted grid results agree 
well, with more small-scale structure in the non-linear stages in the 
meshless methods (owing to increased grid noise). Quantitatively, 
the total magnetic energy in the box grows, with excellent agree¬ 
ment between these methods and converged solutions until t « 2.5 
(well into the non-linear phase). However, the v x = 10 boost com¬ 
pletely suppresses mode growth in stationary grid methods at this 
resolution. In SPH, a reasonable, answer can be obtained with suf¬ 
ficiently high neighbor number and an appropriate choice for the 
artificial resistivity, but the instability is totally suppressed if we 
use typical neighbor numbers and behaves as if the field were much 
weaker in the “vf ast AR” case. Comparing V ■ B, we see the same 
behavior around the contact discontinuity as in the RT test; the rel¬ 
ative performance of different methods is essentially identical, al¬ 
though in all cases the V ■ B errors are systematically smaller by a 
factor of a couple. 

Finally, in Fig. [22] we also compare an MHD version of the 
“blob” test from Agertz et al)] |2007[ . The setup is identical to our 
detailed study of different methods on this problem in Paper I (fea¬ 
turing a cold cloud in pressure equilibrium in a hot wind tunnel), 
but with an initially constant field in the vertical direction. As ex¬ 
pected the field strongly suppresses the non-linear RT and KH in¬ 
stabilities that tend to disrupt the cloud in the hydrodynamic case 
(for a detailed study, see |Shin et al.|20d8| . Our qualitative conclu¬ 
sions (comparing different methods) are identical to the tests above. 

3.10 Magneto-Rotational Instability: Can Meshless Methods 
Capture the MRI? 

We next consider the magneto-rotational instability (MRI), one of 
the most important and historically challenging MHD phenomena. 
We consider a 2D axisymmetric shearing box defined as in |Guan| 
|& Gamm?e| l j2008| ; this is a locally-Cartesian box where the x co¬ 
ordinate represents the radial direction R , and the other coordi¬ 
nate is the vertical direction z (azimuthal axisymmetry is assumed). 
Boundary conditions are periodic in z and shear-periodic in x: 
f(x, z) = f{x+n x L x , z + n z L : ) for all values / except the azimuthal 
velocity v# = v y (x , z) = v y (x + n x L x , z + n z L z ) + n x qQ.L x , where 
q = — (l/2)r/lnfI 2 /<71nR = 3/2 for a Keplerian disk, n x and n y are 
integers representing the box periodicity, and D is the mid-plane 
orbital frequency. In this approximation, the momentum equations 
are also modified with the source terms D(p\)/Dt = —2(f 1z) x 
(p\) + 2pqQ 2 xx. We initialize a box of unit size (—0.5 < x < 0.5, 
—0.5 < z < 0.5) and Q = 1, with (p, P, v x , v y , v z , B x . B y . B z ) = 
(1, 1,0, Sv, —qx, 0, 0, Bq sin(27rx)), where Bq = vT5/(87r m) is 
set so that the most unstable wavelength Amri = 1 /m corresponds 
to a mode number m. Here Sv is set to a uniform random number in 
the range ±0.005 for each cell, to seed the instability. 

shows the results at various times from a 256 2 MFM 


Fig. 


23 


calculation with m = 4. We see the expected behavior; the random 


© 0000 RAS, MNRAS 000, 000-000 






























Meshless MHD 17 


MFM 

MFV 

Grid 

Grid+boost 

SPH-hi 

SPH-lo 

SPH: Vfast AR 

Powell 

LU 

LU 

LU 


LU 


m 

LU 



1 1.5 2 2.5 









Figure 19. Magnetic Rayleigh-Taylor instability (§ |3.9) . We show the density p (as labeled), in both early (top, t — 6) and late-time ( bottom , t = 16) results at 
128 x 256 resolution. MFM, MFV, and grid methods (with no bulk motion) all converge to the correct solution. In the early stages they are nearly identical; 
even in late stages there is little difference (small differences appear at late times in MFV owing to growth of the slightly-larger grid noise at early times). As 
in the pure-hydro case, in non-moving grid methods, “boosting” the fluid grid by a uniform velocity slows down mode growth and breaks symmetry, unless 
significantly higher resolution is used. In SPH, reasonable results can be obtained if large neighbor numbers are used (although some large-scale asymmetry 
appears in the non-linear solution because of imperfect V • B-cleaning); with smaller neighbor number (the same as used in MFM/MFV methods), low-order 
errors in SPH dominate and no mode grows. However, the SPH results are also very sensitive to the form of the artificial dissipation (artificial resistivity or 
AR) employed: we compare the results if we use the fast magnetosonic velocity Vf as t in AR, as advocated by |Tricco & Price](2013) (which reduces noise in 
super-sonic problems), instead of our default choice (the Alfven velocity). For details, see AppendixjA] This produces too much dissipation around the contact 
discontinuity, which damps the B-field to sub-critical values; so the instability behaves (incorrectly) like the pure-hydro case. 



Figure 20. Divergence errors in the MHD RT instability from Fig . 1 1 9| We 
plot log 10 (/z|V • B|/|B|) (as labeled), for the same times as Fig. |19| The 
grid methods here use CT so V • B = 0 to machine precision. MFM/MFV 
methods maintain log 10 (/i|V - B|/|B| ) < —2 (for every particle/cell) even 
into the late non-linear evolution (most of the values at > — 4 are boundary 
effects, in fact). In SPH, the zeroth-order errors around a contact discon¬ 
tinuity lead to less-accurate cleaning in these regions. In the Vf as t artificial 
resistivity (AR) case, /z | V • B|/|B| ~ 1 around these discontinuities; this is 
not because V • B is large, but because the AR over-damps B at the discon¬ 
tinuity. In the Powell case, order-unity errors appear at late times. 

B-field fluctuations seeded by the velocity perturbation grow, and 
quickly are dominated by the fastest-growing (m = 4) mode, until 
at late times the non-linear modes break up into MRI turbulence. 
The same behavior appears in MFV, grid, and even SPH methods. 



Figure 21. Magnetic Kelvin-Helmholtz (KH) instability (§ |3.9| . We show 
density p (as labeled) in the non-linear phases t = 1.6 and t = 3.2 (KH 
growth timescale « 0.7), at 256 2 resolution. Differences between methods 
resemble the RT test (Fig. |T9). MFM, MFV, and non-boosted grid results 
agree well into the non-linear growth phase (where different grid noise ef¬ 
fects lead to small departures). Boosting the fluid in the grid case or using 
typical neighbor numbers for SPH break symmetry and suppress the insta¬ 
bility at this resolution. Divergence errors in each case closely resemble 
those in Fig. |20| around the contact discontinuities. 

In Fig. [24] we compare results when the linear mode dominates, 
for runs with initial B- set such that m— 1, 2, 4, 8, with MFM at 
the same resolution. In each case, we see the correct mode grows. 
Fig. [25] compares different methods and resolutions (for m = 4). 
We see that MFM and MFV produce nearly identical results. Even 
at 32 2 resolution, a reasonable mode structure emerges. In general, 
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Figure 22. Magnetic “blob” test (§ |3.9| . We show density log 10 (p) (arbi¬ 
trary units) at t = 2 and t = 4, with initial vertical field Bq = 10 —5 (top) 
and 10 _1 ( bottom ) in code units. The weak-field case is essentially identi- 
cal to the pure hydrodynamic case studied in Paper I, as expected; a com¬ 
bination of non-linear Rayleigh-Taylor and Kelvin-Helmholtz instabilities 
“shred” the cloud. The strong-field case suppresses cloud destruction. MFV 
is shown, but differences between methods are qualitatively identical to the 
RT and KH tests above. 


with MFM/MFV methods, we find about ~ 4 resolution elements 
(particles/cells) in linear dimension across each “node” are needed 
to see reasonable modes (for m = 2, we see growth for 16 2 reso¬ 
lution); this is the same as in higher-order (PPM) grid-based codes 
using CT (see |Guan & Gammie|2008) . 

For SPH, we see - perhaps surprisingly - reasonable behavior, 
if we use large enough neighbor numbers. Compared to previous 
SPH MHD implementations, the fact that the new methods com¬ 
bine many basic hydrodynamic improvements (larger kernels to re¬ 
duce noise, evolution of conserved variables, Lagrangian-derived 
magnetic forces, and properly-derived divergence-cleaning) leads 
to the ability to follow the MRI. A very similar implementation 
of SPH MHD by Vanaverbeke et ah] ( |2014| > and preliminary re¬ 
sults from |Tricco|j2015| have also demonstrated this. However the 
noise level in SPH is higher, even with this larger neighbor num¬ 
ber, at lower resolutions. And if we use a smaller neighbor number 
in SPH, comparable to what is used in our MFM/MFV methods, 
we see the SPH gradient errors totally dominate the solution, and 
prevent any mode growth (see also |Tricco|2015j l. 

Fig. 26 compares h\ V • B|/|Bj. Even in the late non-linear tur¬ 
bulence stage, these are maintained at < 1CP 2 (median ~ 10~ 4 ) in 
MFM/MFV and < 0.1 (median ~ 1(T 3 ) in SPH. 

Figs. |27|28| compare the growth of the MRI modes quanti¬ 
tatively. We measure the amplitude of the m = 4 Fourier mode 
in radial (B x ) or azimuthal (B y ) field components (the maximum 
of the vertical m = 4 mode amplitude in the 2D Fourier trans¬ 
form), at each time, for the simulations above. We also compare 
the volume-averaged magnetic energy {Eb) = (|B| 2 /2) in the box. 
The linear-theory prediction is that the fastest-growing B x and B y 
modes should grow oc exp(0.75t), so the magnetic energy should 
scale (Eb) = Eo + SE exp (1.5 t), where Eq is the initial energy and 
5E is a seed perturbation amplitude. For B x , B y , and Eb, we see 
good convergence with MFM/MFV to the correct linear growth rate 
at resolutions above ~ 64 2 for m = 4; this agrees well with state-of- 
the-art CT grid methods. At late times, the mode saturates and then 
the energy must decay according to the anti-dynamo theorem; this 
is expected. In SPH, we see no MRI growth unless we go to large 
neighbor numbers; we then do obtain growth, but the convergence 
in growth rates is slower. 

In Fig. [28] we compare the late-time evolution of the magnetic 
energy to grid methods. The linear growth rate and peak |B| ampli¬ 
tude agree well with high-order grid-CT calculations. The late time 
decay rate of the magnetic energy is known to be sensitive to the 
numerical diffusivity of the method (see|Guan & Gammie||2008[ 


and references therein): we therefore compare three different grid 
codes at the same resolution (and one much higher): HAM (most 
diffusive), ZEUS, and ATHENA (in PPM, CT, CTU mode; least dif¬ 
fusive). The late-time behavior in our meshless methods agrees re¬ 
markably well with ATHENA. 

3.11 Collapse of a Magnetized Core: Preserving Symmetry 
& Launching MHD Jets 

Next, we consider collapse of a magnetized proto-stellar core and 
launching of a proto-stellar jet. This is a less quantitative problem; 
however, there are several key qualitative phenomena. A rotating, 
weakly-magnetized, self-gravitating gas sphere is initialized. This 
collapses under self-gravity quasi-spherically to much higher den¬ 
sities, testing the ability to follow the fluid in compressions and 
collapse over many orders of magnitude, and whether the gravity- 
hydro coupling is conservative. The collapse is arrested by the for¬ 
mation of a disk (requiring good angular momentum conservation 
in the code), which slowly contracts via magnetic braking. The 
braking and field amplification in collapse and the subsequent disk- 
driven dynamo test the ability of the code to follow initially weak 
fields and MHD instabilities. Finally, a jet is launched: following 
this requires the code have good symmetry preservation, and most 
critically maintain low values of V ■ B, or else the protostar will be 
ejected from the disk by accumulated errors (see |Price et al.|2012jl. 

Following Hennebelle & Fromang (20081, we initialize a 3D 
box of size-length 0.15 pc, in which the hydrodynamic forces are 
periodic but gravity is not. We initialize a constant-density sphere 
of radius Ro = 0.015 pc and mass lMg, in a uniform, non-moving 
background (filling the box) with factor « 360 lower density. The 
sphere is set in rigid-body rotation such that the orbital time is 
4.7 x 10 5 yr; this corresponds to a ratio of kinetic-to-potential en¬ 
ergy KE/|PE| « 0.045 « 1. The magnetic field is initialized to a 
constant value Bq, aligned with the angular momentum vector of 
the sphere. At all times, the system is forced to obey the barotropic 
equation of state P = (0.2kms _1 ) 2 p\/l + (p/po) 4 / 3 , with po = 
10 _14 gcm _3 . The value of the field Bq is chosen to correspond 
to different mass-to-magnetic flux ratios, defined in the traditional 
fashion relative to the “critical” value quasi-stability of a spherical 
cloud, p = (M/<E>)/(M/4?) cr i t where (M/<5>) = M/(nRlBo) and 
(in our units) (M/$) cr it = 0.126G -1 / 2 . With M fixed, p oc B^ 1 , 
and for our choices p = 10 corresponds to Bo = 61 pG. 

Fig. [29] shows the results from our MFM method, for various 
values of p (So), using 60 3 total cells/particles in the box, or « 50 3 
cells in the collapsing sphere (we use equal-mass particles, so the 
initial packing is denser in the sphere). The times are chosen to 
be some (short) time after a jet forms, in each case close to t ~ 
4 x 10 4 yr ~ ts, where ts = %/3/(2-izGp) is the gravitational free- 
fall time. Fig.[3T]shows the same for MFV. 

Even at this (modest) resolution, we are able to see all of 
the key behaviors above over a wide range of p. As expected, 
for p < 1, a stable, thick disk forms, which does not spin up the 
field sufficiently to launch a jet. But for larger p > 1 (where the 
disk collapses rapidly), a jet is launched and propagates well into 
the diffuse medium as accretion onto the protostar continues. The 
collimation and strength of the jets vary with p as expected, and 
are consistent with AMR simulations using CT and run with much 
higher resolution (see jHennebelle & Fromang|2008| >. A more quan¬ 
titative exploration of this is an interesting scientific question, but 
outside the scope of our study here. We stress that as the fields be¬ 
come weaker, this problem becomes more challenging, since any 
excess numerical dissipation will tend to suppress field amplifica¬ 
tion and jet-launching. This is well-known from grid-based studies 


© 0000 RAS, MNRAS 000, 000-000 

































Meshless MHD 19 



Figure 23. Magneto-rotational instability (MRI; § |3.10) test problem. We show the results from a 2D axisymmetric shearing box (horizontal/vertical axes 
correspond to radial/vertical coordinates), with vertical B z set so the fastest-growing modenumber is m = 4. We plot the azimuthal/toroidal component B y of 
the magnetic field (the behavior of the radial component is similar, but with reversed sign), in an MFM calculation at 256 2 resolution (units are scaled to the 
maximum/minimum values in each frame as labeled, since the absolute value of B y grows exponentially). The initial seed noise is amplified on the correct 
timescale, and from times t ~ 5 — 16, the dominant m = 4 mode pattern is clearly visible. At late times, the non-linear modes break up into turbulence, as 
expected. MFM, MFV and high-order CT-based grid methods (see |Guan & Gammie|2008) produce similar results. 



Figure 24. MRI, as Fig. |23| (same colorscale; just showing the 0 < x < 
0.5 half of the box), at times when the linear mode dominates, for runs 
with initial B z set so that the fastest-growing mode corresponds to m = 
1.2,4, 8, as labeled (MFM shown, MFV is nearly identical). In each case, 
the correct mode is clearly visible. We find approximately 4 cells/particles 
are required across each “node” in the linear direction to resolve the correct 
mode growth, the same as the number of cells in PPM grid methods. 

using Riemann solvers with different diffusivity. What is remark¬ 
able is that we are able to see jets even at /r > 10 at this resolution. 
For weaker fields than about /j, ~ 20, various authors have noted 
that the magnetic field becomes insufficient to stabilize the disk, 
which then fragments; we see this in our weak-field case. Flowever, 
even here, we see that each fragment/protostar is launching its own 
“mini-jet,” which then precess owing to the orbital motion of the 
fragments! 

To see these behaviors, particularly in the weak-field case, 
usually requires quite high resolution (> 128 3 ) in grid-based codes 
i jHennebelle & Fromang|2008]|Pakmor et al.|201 1| >. In Fig. [30] we 
show a resolution study at fixed time in our fiducial MFM, /r = 10 
case. While the disk continues to show more fine structure, and 
the jets are launched slightly earlier (so have propagated further) 
at higher resolution, the presence of a jet is clearly resolved, and 



Figure 25. MRI (with mode number m = 4), as Fig. [23] as a function of 
resolution and numerical method. MFM/MFV methods are very similar, 
and resolve the MRI with m = 4 with as few as 32 2 elements. Modern SPH 
MHD is, in fact, also able to capture the MRI modes, if a larger neighbor 
number is used. If we use a smaller neighbor number in SPH, the low-order 
errors completely swamp the correct solution. 

the jet outflow rate and momentum are converged to within tens of 
percent, at 25 3 resolution; some jet is even visible at 12 3 resolution! 

We can see these behaviors even at very low resolution be¬ 
cause this test combines many advantages of our new methods. It is 
clearly Lagrangian, with high-dynamic range collapse; the fact that 
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Figure 26. Divergence errors in the MRI test in Fig. |25| at 256 2 resolution. 
We plot log 10 (/t|V • B|/|B|) (as labeled), for two times, t = 10 when the 
fastest-growing mode dominates, and t = 20 when the system has broken 
up into turbulence (Fig. [23). In both cases, the errors are well-controlled 
in MFM/MFV; in SPFI some regions reach larger h | V • B|/|B| > 0.01 in 
the turbulent phase; however these almost entirely correspond to regions of 
nearly-vanishing |B| (i.e. (h | V ■ B|)/(|B|) <7 0.01), so the errors are not 
dynamically significant. 


(V-body codes naturally couple to our methods with fully adaptive 
gravitational softening providing second-order accuracy (Paper I) 
means that we can follow the entire collapse without any special 
“switches” required in our gravity solver. The smooth, continuous 
adaptive resolution provided by our method avoids the low-order 
errors which tend to artificially damp magnetic fields, inherent to 
AMR refinement boundaries. Once the disk forms, it is critical that 
the method be able to accurately conserve angular momentum and 
not degrade orbits; this is especially challenging in AMR codes (as¬ 
suming ~ 50— 100 3 cells across the central disk, even high-order 
AMR methods will degrade the orbits within a couple of periods; 
see Paper I &|de Val-Borro et al,|2006[|Hahn et al.|20T0]|Duffell &| 
|MacFadyen|2012) . 

Fig. 32 demonstrates that, because they are Lagrangian and 
mesh-free, our methods here are trivially invariant to both (a) boost¬ 
ing the entire fluid (so the core is moving super-sonically through 
the box), and/or (b) rotating the core at an arbitrary angle to the co¬ 
ordinate axes. This is not true in non-moving grid methods: either 
of these changes in grid codes implies a very significant loss of “ef¬ 
fective resolution” or accuracy at fixed resolution. Based on com¬ 
parison experiments with RAMSES, we estimate that for the p = 10, 
rotated+boosted case in Fig. [32] achieving comparable accuracy in 
AMR to our ~ 50 3 (10 5 -element) MFM run requires ~ 10 7 8 cells. 

SPH, however, has greater difficulty with this problem. Fig.|33| 
shows the same survey of p at 50 3 resolution in SPH; in no case 
does a stable jet develop. For small initial fields (p > 10), artifi¬ 
cial resistivity too-efficiently dissipates the field, so the case resem¬ 
bles pure hydro (either weak or no jets appear). For large fields, 
the method has difficulty maintaining stability: the central disk 
breaks up, and the proto-star is ejected! These errors are resolution- 
dependent, however. At higher resolution, we do begin to recover 
the correct behavior in SPH; Fig. 


34 


shows the same at 10(T reso¬ 
lution. Here the high-B cases begin to resemble the MFM/MFV re¬ 
sults. The low-B cases still exhibit too much dissipation (too-weak 
jets), but the field values are (very slowly) converging towards the 
MFM/MFV resultE] 


19 We have explicitly checked, for the proto-stellar jet test in SPH, whether 
using the Alfven or magnetosonic velocity in the artificial resistivity (as 


Fig. [35] compares the divergence errors; this illustrates the 
source of the instability in SPH. In MFM/MFV (essentially iden¬ 
tical here), the errors are well-controlled at the 10 -3 — 10~ 2 level 
(a few particles in the disk midplane reach larger values, but these 
are in the regime where the thin disk is completely unresolved ver¬ 
tically so the change in B must occur across a single particle). We 
also clearly see that the errors decrease with resolution (the me¬ 
dian h | V ■ B|/|B| decreases by a factor of ~ 5 from 50 3 to 100 3 
resolution). In SPH, we see much larger V ■ B, as in many cases 
above - most critically, we see that where the disks have gone un¬ 
stable and “kicked out” the protostars corresponds (in every case) 
to /i| V -B|/|B| 3> 1. This is consistent with previous studies (Biir- 


|zle et a l. 2011; P rice et al.|2012| >. Because the V • B errors decrease 
with resolution, going to higher-resolution suppresses this and al 
lows stable evolution of the jets. 


3.12 MHD Zeldovich Pancake: Testing Cosmological MHD 
Integration & Extremely High Mach Number Shocks 

We now consider an MHD Zeldovich pancake, following a sin¬ 
gle Fourier mode density perturbation in an Einstein-de Sitter 
space, to test the code’s ability to deal with cosmological inte¬ 
grations, small-amplitude perturbations, large Mach numbers, and 
anisotropic cell arrangements. We initialize a linear perturbation 
following jZeTdovichj |1970|, with parameters chosen to match 
|Li et al.| (2008$ , in a periodic box with side-length L = 1 and 

7 = 5/3. Assume the unperturbed fluid elements have uniform 

density with co-moving position q along x as z —> oo, then co¬ 

moving positions and fluid quantities at the initial redshift z, = 20 

of the simulation are* = q — As'm(kq)/k,p = po/(l—A cos (kq)), 

Vpec = —H 0 sin [kq)/k, where A = (1 +z c )/(l +zi),k = lit/L, p 0 = 
3Ho/(SnG) the critical density, and z c = 1 is the redshift of caustic 

formation. In our particular (arbitrary) code units. Ho « 0.18, and 

the initial u = 9.3 x HP 8 and B = Boy (comoving). We consider a 
weak and strong field case, with B 0 = 1.25 x 10 -4 and 7.7 x 1CP 3 , 
respectively. These are rather unconventional units, but we note 
that the solution can trivially be rescaled; it depends only on the 
dimensionless parameters z c and the initial /3 = Rthermai//magnetic- 
The perturbation should grow linearly along the x direction under 
self-gravity, until it goes non-linear and eventually collapses into a 
shock (caustic) at z = 1. 

At early times, the flow is smooth and obeys a known linear 
solution; we confirm that our MFM/MFV methods reproduce this 
with second-order convergence (even in cosmological integration 
with self-gravity coupled to MHD). In Fig.[36]we plot the results at 
z = 0. There are now clear strong shocks (pressure jumps of factors 
~ 10 s ) and factor ~ 1000 compressions. At various locations, the 
kinetic, thermal, and magnetic energy form an extremely disparate 
hierarchy: in the weak-field case p|v ra m| 2 ~ 10 lo * * * pc 2 ~ 10 14 |B I 2 , or 
in the strong-field case, p |v ra m| 2 ~ 10 9 |B| 2 ~ 10 15 pc 2 . As a result, 
Eulerian codes which evolve only the total energy almost invariably 
crash or produce unphysical results (e.g. subtraction of very large 
numbers giving negative pressures). The Lagrangian nature of our 


in § |3.9} makes a large difference; it does not. Likewise the choice of 
“pressure”-SPH (PSPH) versus traditional “density”-SPH (TSPH), dis¬ 
cussed in Paper I, is not especially significant for the core/jet problem. We 
have also experimented with a wide range of different artificial viscosity and 
conductivity parameters and find the qualitative results are similar. How¬ 
ever, as noted in |Tricco & Price] (2013) , it is important that the artificial 
resistivity be limited to a small maximum value; our default here is 0.1 (see 
Appendix[Aj, if we raise this to 1, the B fields are artificially damped away 
to negligible values during collapse. 
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Figure 27. Growth of the m = 4 MRI modes in Figs. |23|25| We measure the m = 4 Fourier amplitude in the azimuthal (B y ; left) and radial (B x \ center) 
magnetic fields, as well as the average magnetic energy ((|B~/8vr|) in the box (right). We compare each to the expectations of linear theory (dotted black line). 
We compare three different resolutions (64 2 , 128 2 , 256 2 , in progressively thicker lines). MFM/MFV methods give similar results; in both cases the simulations 
converge to the correct linear growth rate quickly at resolutions above ~ 64 2 . This agrees well with state-of-the-art grid CT codes. The peak mode amplitude 
also increases with resolution; once the MRI breaks into turbulence, the field amplitude declines, as it should. For SPH, we see no MRI at any resolution unless 
we use a large neighbor number. With sufficiently large neighbor number, we obtain growth, but convergence to the correct growth rate is slower (even at 256 2 , 
the growth rate is suppressed by ~ 30%). 



Figure 28. Late-time evolution of the magnetic energy, for the m = 4 case, 
at 256 2 resolution. We compare results from three different unsplit, CT- 
based grid codes in |Guan & Oammlc|(20()8} . The 2nd-order HAM code is 
most diffusive, while the 3rd-order PPM ATHENA is least-diffusive (ZEUS 
is intermediate). In all cases, the linear growth rates and peak amplitudes 
are similar and agree well with our meshless methods; however the decay 
rate is sensitive to the details of numerical dissipation. Our meshless results 
at fixed resolution are most similar to ATHENA, the least-diffusive of the 
grid-based codes considered. 

method greatly reduces these errors, and the dual-energy formal¬ 
ism (Paper I; Appendix D) allow the method to deal smoothly with 
disparate hierarchies. 

Another key aspect of the 3D test is that it involves highly 
anisotropic compressions: the gas collapses by a factor of ~ 1000 
along x, but there is no collapse in the other directions. We discuss 
this at length in Paper I: for AMR methods, this makes it very ex¬ 
pensive to achieve the same resolution as our mesh-free methods, 
because (since grid cells are cubical), refinement must take place 
along the y and z dimensions with the collapse; for moving-mesh 
methods, it requires very careful and accurate cell refinement and 
regularization schemes (which can introduce other errors). But it 
is handled continuously and naturally in our mesh-free Lagrangian 


schemes. For the same total number of resolution elements, this al¬ 
lows our Lagrangian methods to achieve a factor ~ 10 better spatial 
resolution along x in the central regions, compared to AMR (in 3D). 

In the weak-field case, we confirm the results from Paper I: 
the method can handle large energy hierarchies; arbitrarily strong 
shocks are resolved across ~ 2 linear cells/particles for MFM/MFV, 
and smeared across a factor ~ 2 — 3 larger range in SPH; co-moving 
integration with self-gravity is accurate; and because of their abil¬ 
ity to handle asymmetric cell/particle distributions, the Lagrangian 
methods converge more rapidly in high-density regions compared 
to AMR. In the strong-field case, we also confirm that trace B-fields 
are correctly amplified, and recover the correct jumps and rarefac¬ 
tions. 

On this particular test, we obtain similar results with the |Pow-| 
|ell et al.| ( |1999[ l-only divergence cleaning. The reason is that there is 
negligible field in the x direction, and the growth of the perpendic¬ 
ular By component is driven by simple, effectively one-dimensional 
compression/expansion. 

3.13 The MHD Santa Barbara Cluster: Cosmological MHD 
Integration in Turbulent Flows & Divergence-Control 

Next we consider the MHD “Santa Barbara Cluster” from {Frenkj 
|et al.|[T999] >; the hydrodynamic case is again in Paper I. The test 
is a “zoom-in” where we initialize a high-resolution Lagrangian 
region in a low-resolution Einstein-de Sitter cosmological back¬ 
ground, which collapses to form an object of a rich galaxy cluster 
mass at z = 0. The cluster ICs are in jFrenk et al.| < fl999| l; a peri¬ 
odic box of side-length 64 h~' Mpc is initialized at redshift z = 49, 
in a flat Universe with dark matter density Hdm = 0.9, baryonic 
S2b = 0.1, Hubble constant Ho = 50 km s - 1 Mpc -1 . The gas is non- 
radiative with 7 = 5/3 and initial T = 100 K. We add to this a trace 
initial seed field B = Bqz- The initial B direction and magnitude 
should be unimportant, provided it is small. 

Fig. [37] shows the resulting profiles of density, temperature, 
magnetic field, and V ■ B at z = 0, across simulations using different 
methods (the average jB| and h \ V ■ Bj/|B| values plotted in Fig.|37| 
are magnetic energy weighted averages). The grid result here is 
taken from high-order (PPM unsplit CTU) AMR simulations using 
CT in |Miniati & Martin| l |2011| >. The dark matter density profile (es¬ 
sentially determined by the (V-body solver) is nearly identical in all 
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Figure 29. Proto-stellar jet test (§ |3.1l) . A magnetized, rotating proto-stellar core with initial mass-to-flux ratio p oc |B| 1 collapses under self-gravity until 
it forms an accretion disk, which winds up the magnetic field and launches a jet. We show results for MFM calculations with 50 3 total particles/cells in the 
core, and various p (columns labeled), just after one free-fall time of the initial core. Scale bar in each panel shows 200 au. Top: Density (log 10 («/cm _3 )), 
as labeled) in a slice through the collapsed disk midplane (face-on). The central “protostar” has collapsed by a factor of > 10 4 in density; the disk is more 
extended and lower-density for stronger B. For p < 1, the disk is quasi-stable, and should not form a jet. For p > 20, the disk should fragment into multiple 
protostars. Middle: Density (log 10 (/7/cm -3 )) in a slice through the rotation axis. Here the jets are plainly visible. Bottom: Radial velocity (v r /kms -1 ) in a 
slice through the rotation axis. Even at low resolution, and very weak B, our new meshless methods are able to capture jet-launching. 



Figure 30. Resolution study of the collapsing core as in Fig. |29| using 
MFM (MFV is similar), with initial p = 10. We compare at fixed time 
(t = 1.05%). The resolution quoted is the total number of cells/particles in 
the collapsing spherical core. At higher resolution, additional fine-structure 
in the disk and outflow continue to appear; the jet also forms slightly earlier, 
so it has propagated further. However, there is already some (albeit weak) 
jet/polar outflow structure at 12 3 ; by ~ 25 3 , the existence of a jet and much 
of the global structure is already resolved, and the jet momentum/mass are 
converged within tens of percent of the highest-resolution run. This is re¬ 
markably low resolution; usually, > 128 3 resolution is needed for similar 
convergence in grid calculations. 

runs, as expected (as is the final pressure profile, given by hydro¬ 
static equilibrium). There are some very small differences in the 
central gas density, temperature, and entropy; these are discussed 
extensively in Paper I. 

In all methods, for a wide range of Bo, B is amplified to ~ pG 
in the cluster core. In our “default” runs, we seed Bo = 10 -8 G 
(4 x 10“ 12 G co-moving). This gives plasma p = Pgas/Pmag > 50 
everywhere in the ICs, so the magnetic pressure is unimportant. 
As long as this is true, the final B profile is nearly independent 
of Bo . We have verified this in all methods: for Bo > 10 -7 G, the 



results are similar up to non-linear details. Note that the p = 20 case is 
marginally unstable to fragmentation; so small changes to the particle split¬ 
ting/merging algorithm in MFV (which can seed grid noise) can lead to 
modest perturbations that induce earlier fragmentation (and weaker jets); 
the “default” case here allows no splitting/merging (most similar to MFM). 


initial ft ~ 1 and the hydrodynamic properties (p, T) as well as 
maximum |B| begin to change. For sufficiently small Bo, in prac¬ 
tice, numerical resistivity and truncation errors will swamp the field 
and suppress growth: we obtain similar final B profiles for mini- 
mum Bo > (10~ 14 , 10 -11 , 1(T 9 , 10“ 12 ) G in (MFM, MFV, SPH, 
AMR), respectively. The fact that we can use such small minimum 
Bo (comparable to roundoff errors) in MFM reflects the extremely 
low numerical dissipation of advected quantities inherent to the 
method. The large value in SPH reflects the artificial resistivity er¬ 
rors discussed above, whereby fields with /3 100 tend to be arti¬ 

ficially over-damped. Even for Bo = 10~ 8 G, we already see some 
SPH over-damping in the cluster outskirts. 

In MFM/MFV methods, V • B is well-controlled, with mean 
absolute values of h | V ■ B|/|B| ss 0.01 in the resolved region of the 
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Figure 32. Effects of rotations and boosts on our meshless methods in 
the protostellar jet test. Left: Our “default” MFM run from Fig. [29] with 
p = 10, resolution 50 3 , zero net bulk motion and an initial angular mo¬ 
mentum/magnetic field axis aligned with the z direction. Right: Same, but 
with the rotation/field axis rotated 30° in the x — z plane, and the entire 
box boosted by a uniform 8v x = lOkms -1 , 5v-y = 5v z = 2kms -1 . Our 
mesh-free methods are invariant to boosts (bulk motion) of the box and ro¬ 
tations/grid mis-alignments at all resolutions; however these pose serious 
challenges for grid-based calculations. 



Figure 33. As Fig. [29] for SPH-hi. SPH has much greater difficulty captur¬ 
ing the important behaviors at low resolution. At low-|B|, the SPH “artifi¬ 
cial resistivity” over-damps B and suppresses jets. At high-|B| (p <C 10), 
less-accurate V • B-cleaning leads to unstable errors that disrupt the disk. 



Figure 34. As Fig.[33] for SPH-hi but at higher resolution (100 3 ). At higher 
resolution, SPH is better able to control V ■ B errors and jets emerge at 
high-B {fi < 10); however the SPH numerical resistivity still over-damps 
the weak-field case. 



show face-on (upper) and edge-on (lower) slices as Figs. |29|34| for MFM 
(fop; MFV is similar) and SPH-hi (lower panel). The unstable behavior in 
SPH at low resolution is clearly related to poor control of the V • B-errors. 
We show two p values at 50 3 resolution and compare // = 5 at higher(100 3 ) 
resolution. In all cases V ■ B errors decrease with resolution. 


cluster. The errors rise towards the outskirts owing to (1) decreasing 
resolution there, and (2) boundary effects (the setup of this “zoom 
in” involves no gas outside the initial Lagrangian region, so there 
are vacuum boundaries on one side of many particles). The median 
h | V ■ B|/|B| is ~ 10 -4 — 10 -3 . In SPH, the errors are less well- 
controlled, reaching ~ 0.1; however this appears to have little or 
no effect on the solution. 

With Powell-only cleaning, however, the field is artificially 
amplified to order-of-magnitude larger values in the cluster core; 
this is essentially the compounded version of the erroneous growth 
seen in the field loop test. This is qualitatively different from 
the behavior seen in any runs with cleaning, despite the fact that 
A|V-B|/|B| (while large) is not much larger than our with-cleaning 
SPH runs. Once again, this emphasizes that in Powell-only clean¬ 
ing, the non-linear error terms can integrate unstably (i.e. build co 
herently). In contrast, with a properly-applied [Dedner et al.|(2002 


divergence-damping, these terms are stabilized and Tricco & Price] 
( |2012| > show that the B field cannot artificially self-amplify (rather, 
the sense of the errors will be to produce additional numerical dis¬ 
sipation, just like most other sources of error). So even if V ■ B 
is nominally the same over the course of a simulation, runs with 
|Dedner et al. 1(2002) cleaning will avoid many of the more serious 
instabilities of Powell-only cleaning. 
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Figure 36. Magnetic Zeldovich pancake (§ |3.12) at z = 0. An initially small (linear) sinusoidal density perturbation collapses along one dimension (xj in a 
3D cosmological expansion, forming a caustic at z = 1, with a small initial perpendicular magnetic field B = Boy. We plot velocity in the collapse direction 
(i’ A ), density p (relative to the box-mean, po), pressure P. and field B y . All cells/particles are plotted; the resolution is ~ 100 3 . Left: Weak-initial field case. 
Here the results are similar to the pure-hydro case (see Fig. 30) in Paper I; B y grows via compression in the center and declines adiabatically with the Hubble 
flow outside the shock. Right: Strong-initial field case. Here, the shock center is magnetically dominated, flattening the density peak and generating a pressure 
cavity. In both cases, MFM/MFV converge well to the exact solution (here, a 2048 3 -equivalent, one-dimensional grid-based PPM. CT calculation from ENZO) 
- a ID, 512-element MFM calculation is indistinguishable from the exact solution. At lower (100 3 ) resolution, the outer, lower-density shock at x = ±0.03 is 
broadened by ~ 2 particles in the linear direction (about the same as obtained in second-order grid methods). The inner region, including the detailed shape 
of the central density and B-field, are already well-converged. The shape and boundaries of the low pressure cavity in the strong-field case are well-resolved, 
but the lowest value of the thermal pressure converges slowly in the central region where ^thermal "C |B | 2 /2 (i.e. where it is dynamically insignificant). In 
SPH, shocks are broadened by about twice as many particles in the linear direction, and the density peak shape converges slightly more slowly. The shape of 
the boundaries of the central pressure cavity in the strong-field case converge significantly more slowly than in MFM/MFV/grid codes. The V • B errors are 
negligible in this problem. 


3.14 Simulated Galaxies: Testing Code Robustness in 
Multi-Physics Applications 

We now consider a “stress test” of the methods here, adding mag¬ 
netic fields to simulations of galaxies using state-of-the-art physics 
models. Specifically, we consider three initial conditions: (1) an 
isolated (non-cosmological) galaxy (with a pre-existing stellar disk 
and bulge, gas disk, and dark matter halo), designed to represent 
a starburst/M82-like system (model Sbc in |Hopkins et al.|2012c[ 
with a IfiG seed field); (2) a cosmological “zoom-in” simula- 


der to maximize the numerical challenge. For each galaxy, we ac¬ 
tivate the full suite of physics from the FIRE (Feedback In Real¬ 
istic Environments) simulations described in Hopkins et al.|(20T~i~| 


|2012b| |2013a| > and |Faucher-Giguere et al.| ( |2015| >. This includes: 
self-gravity; gas cooling and chemistry; star formation; cosmologi¬ 
cal expansion; the interaction of gas, stars, and dark matter; energy, 
mass, momentum, and metal injection from supernovae and stel¬ 
lar winds; and radiation-matter interactions in the form of photo¬ 
ionization, photo-electric heating, and radiation pressure. 

We emphasize that there is no simple “correct” solution for 


a dwarf galaxy (with z = 0 halo mass 10 l0 MfVi; model 

the B-field evolution in these tests, and our concern here is not 

Hopkins et al. 2014 with a 10~ 10 /rG seed field); and (3) 

whether this particular model of the physics is correct or complete 


mlO in 

a “zoom-in'j_of_a Milky Way-like system, run to z = 2 (model 
ml2i in Hopkins et ahj2014; \0~ >o p,G seed field). We intention¬ 


ally study low-resolution versions of each - using a factor of ~ 64 
fewer particles than the “production runs” in those papers - in or- 


(we know, in fact, that these examples are under-resolved). Rather, 
we test (1) whether or not the algorithms we have developed can 
run (at all), without crashing or returning unphysical solutions; and 
(2) how well they control the divergence errors. This is extremely 
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Figure 37. Magnetic Santa Barbara cluster (§ |3.13) . A high-resolution sub¬ 
volume of a cosmological simulation is followed; it forms a cluster-mass 
dark matter halo; we show radially averaged profiles at z = 0. Top Left: Gas 
density (symbols) and dark matter density (lines). Bottom Left: Tempera¬ 
ture. SPH shows slightly lower central-7" (see Paper I). Top Right: Magnetic 
field. The central rms |B| ~ pG. independent of the numerical method. With 
Powell-only cleaning, however, the magnetic energy is artificially amplified 
(as in the field loop problem) to order-of-magnitude excessive values. In 
SPH, excess artificial resistivity leads to some excessive damping at large 
radii. The “low seed" run features a seed field a factor ~ 10 5 weaker than 
the default run: the final B is nearly independent of the seed. Bottom Right: 
Divergence errors. Absolute values in | V • B| are taken before averaging: 
the mean |(/tV • B/|B|)| ~ 10 s is approximately 6 orders of magnitude 
smaller. In MFM/MFV, the well-resolved (< 1 Mpc) region maintains small 
divergence errors. In SPH the errors are factor ~ 5 larger. 

challenging: essentially every numerically difficult situation our 
test problems have considered above will occur here (and often be 
poorly-resolved). In addition, gas is dis-continuously added & re¬ 
moved from the system (by stellar mass loss and star formation), 
and non-MHD forces (gravity and stellar feedback) are constantly 
re-arranging the particles/cells on timescales much faster than the 
local magnetosonic crossing times. 

Fig.[38]shows the results for our MFM runs, and Fig. [^com¬ 
pares MFV and SPH for a subset of the initial conditions. With 
the “full” FIRE physics, all the runs develop a multi-phase, super¬ 
sonically turbulent medium, with strong galactic outflows (for de¬ 
tails, see |Hopkins et al.|2012a||2013b|c| >. The B fields are amplified 
to values in very rough equipartition with the thermal+turbulent 
energy of the disk: in cold molecular clouds (T = 100 K, v tur b ~ 
lOkms -1 , n ~ 10— 10 3 cm -3 ), this implies |B| ~ 10 — 100 jiG 
and small plasma ft ~ 10 -2 , while in SNe-heated bubbles with 
T ~ 10 7 K and n ~ 0.01 we find ft as large as ~ 10 — 100. The 
algorithms are stable under arbitrarily long integration. 

In general, we find that V ■ B is reasonable well-controlled, 
with mean values of h\ V ■ B|/|B| ~ 0.03 — 0.1 in MFM/MFV, and 
somewhat larger ~ 0.1 — 0.2 in SPH. Although still within this 
range, the errors are clearly larger in our cosmological zoom-in 
runs. This owes to two facts: they are less well-resolved, and they 
are less dynamically relaxed systems (being perturbed by mergers, 
accretion, etc). However, we stress that these relatively high val¬ 
ues of h |V ■ B|/|B| are totally dominated by the regions where the 
magnetic fields are dynamically irrelevant (i.e. where jBj is small). 
In these regions, the fields are passive, and our divergence cleaning 
scheme has essentially no time to “respond” to the constant, super¬ 
sonic, turbulent particle rearrangement. We therefore instead plot 
/z|v ■ B|/i/|B| 2 + 87rP,hermai = h |V ■ B|/|B| s/l + ft, which com¬ 
pares the V ■ B error to the total MHD pressure (the relevant term 



Figure 38. Galactic disk test (§ |3.14) . We simulate a galaxy disk (gas, stars, 
and dark matter) with full self-gravity, star formation, cooling & gas chem¬ 
istry, and stellar feedback/mass return in SNe, stellar winds, and radiation 
(photo-heating and radiation pressure). We consider four galaxy models: an 
isolated (non-cosmological) starburst disk, with either the sub-grid|Springel| 
& Hemquist 1 2003 1 “effective equation of state” model for the ISM (left: 
this does not explicitly treat the small-scale ISM turbulence and multiphase 
structure), or the “full” FIRE physics modules from |Hopkins et aT~j{20l4j 
(secondfrom left', these explicitly treat the multi-phase, supersonically tur¬ 
bulent ISM), as well as a cosmological zoom-in simulation of a Milky- 
Way mass galaxy run to z = 2 with the full FIRE physics ( second from 
right ; undergoing a major merger at this time), and a cosmological zoom- 
in of a dwarf with the FIRE physics (right: at z = 0). Top: Gas density 
log 10 (n/cm —3 ), in a slice through the galaxy midplane, as labeled. Second 
Row: Plasma log 10 (/J) (ft = Pthermal/Tmagnetick in the same slice. The B 
values reach approximate equipartition with the thermal+turbulent energy, 
producing a wide range in ft. Third Row: Divergence error log 10 (/r|V • 
B| /1 B |). These are reasonably controlled but can reach local values > 
10%, because particles are being re-arranged by gravity and feedback on 
timescales much faster than the fast magnetosonic “response time” for 
divergence-cleaning. Bottom: log 10 (/f | V - B|/-v/|Bp + 8 tt/ jhcrmal )■ This 
shows the divergence error relative to the total hydrodynamic pressure: here 
typical values are < 10 —2 - the large nominal values of | V • B|/|B| gener¬ 
ally appear only in regions where the magnetic fields are dynamically irrel¬ 
evant. 

for the MHD forces); here we see that the errors are actually well- 
controlled at the < 10~ 2 level. It is also worth noting that, like 
in |Pakmor et al.|f2011[ l, the errors are “locally offsetting” - if we 
smooth/average the fields over a few neighboring particles, | V ■ B| 
rapidly vanishes. 

We also compare a more simplified ISM/star formation model; 
this is the sub-grid Springel & Hemquist (20031 “effective equation 
of state” model. Here, the turbulence and phase structure of the ISM 
is not resolved, but replaced with a prescription which forms stars, 
kicks gas out of the galaxy in a “wind,” and pressurizes the gas 
such that certain large-scale properties can be recovered (this is the 
model used in popular large-volume simulations such as Illustris or 
EAGLE; |Vogelsberger et al.|2013fr . By design, in these “sub-grid” 
models the small-scale phase structure and turbulence are smoothed 
over in the ISM (although the simulation still includes super-sonic 
motion, self-gravity, star formation, and galactic winds). With a 
smoother gas distribution, our divergence-control scheme does an 
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Figure 39. As Fig.[38] but for our MFV and SPH-hi methods. Here we just 
compare the isolated disk, with the simplified sub-grid (“smoothed ISM”) 
physics or “full ISM” (multiphase, turbulent) physics, as labeled. The gas 
density distributions, star formation rates, and galactic outflow rates are 
similar (for the same physics) in each method, although SPH predicts some¬ 
what stronger/weaker fields in the inner/outer regions of the galaxy. SPH 
has larger divergence errors by a factor of ~ 3 — 10. consistent with our 
previous tests. 

excellent job maintaining errors < 1(V 2 . Compare this to Figs IS¬ 
IS in |Pakmor & Springel|2012] who apply the Powell-only scheme 
to essentially the same problem in AREPO, and find mean errors 
h | V ■ B | / |B| « 1 (nearly independent of their resolution). 


about 15%. Compared to many other MHD implementations - 
particularly constrained-transport methods - this is extremely ef¬ 
ficient: in ATHENA using CT the speed difference is a factor as 2.3 
(and in moving-mesh methods can be much larger; see |Mocz et al.| 
|2014a} . The memory requirements of MFM/MFV/SPH are also 
somewhat higher in MHD compared to pure hydro, but only by 
the amount needed to carry the additional MHD quantities (B, i/>), 
their gradients, and time derivatives. 

Comparing performance across different MHD methods in our 
tests is non-trivial, because the non-linear dynamics are slightly dif¬ 
ferent. But on all our tests, MFV and MFM are nearly identical in 
cost (MFV is systematically « 5% slower, owing to the need to cal¬ 
culate mass fluxes). If we consider an idealized, pure-MHD test and 
force identical timesteps (i.e. compare cycles per second), we find 
that MFM is only as 5 — 10% more expensive than SPH-lo (which 
requires no reconstruction or Riemann problem and uses the same 
neighbor number); however as in the hydro case, larger timesteps 
are allowed in MFM, so we actually find in our real runs that MFM 
is slightly faster than even SPH-lo and “bare-bones” SPH MHD 
with the same neighbor number. SPH-hi (with IVngb = 120 in 3D), 
on the other hand, is a factor as 2 — 2.5 more expensive than MFM 
(because ~ 4x as many neighbors are needed). 

Even comparing complicated, highly non-linear problems 
with gravity, like the MHD Zeldovich pancake or Santa Barbara 
cluster, we find similar results. On these problems MFV is as 10% 
more expensive than MFM ( the additional difference owes to vari¬ 
able particle masses in the gravity tree), and SPH-hi is a factor 
as 2 — 2.5 slower (for the same particle number). Furthermore, 
since we have shown that achieving the same accuracy in SPH 
requires significantly larger number of resolution elements (and 
convergence in SPH is slow), we conclude that, at fixed accuracy, 
our MFM/MFV methods are less expensive than SPH by factors of 
~ 2 — 10 (depending strongly on the problem). 


4 PERFORMANCE 

In Paper I, we compare performance in terms of speed and mem¬ 
ory useage, across MFM/MFV, our modern SPH-hi, and “bare- 
bones” SPH (the fastest but least accurate form of SPH, using 
no higher-order switches for diffusion, the simplest SPH forms of 
the hydrodynamic equations, and small neighbor numbers), and 
both moving-mesh and grid/AMR codes. Although performance 
is always problem-dependent, we found in general that speed (at 
fixed number of resolution elements) in MFM/MFV was compara¬ 
ble to (or slightly faster than) “bare-bones” SPH, and a factor of 
~ 1.5 — 2.5 faster than SPH-hi, which is itself comparable to the 
speed of moving-mesh algorithms. The memory requirements are 
very similar across MFM/MFV/SPH-hi/SPH-lo, and substantially 
(factor ~ 2) lower than in moving-mesh or AMR methods. 

In terms of zone-cycles per second, the MHD versions of 
the MFM/MFV algorithms are slower than the hydro-only algo¬ 
rithm by about as 30% (run time a factor of 1.3 larger), owing 
to the additional fluid quantities (and their gradients) which must 
be evolved and reconstructed (and extra terms in the Riemann 
solver/equation of motion)]^] The fractional difference in SPH is 


5 DISCUSSION 

We have extended the mesh-free MFM and MFV finite-volume, ar¬ 
bitrary Lagrangian-Eulerian Godunov-type hydrodynamics meth¬ 
ods from Paper I to include ideal MHD. We have implemented a 
second-order accurate, conservative formulation of these methods 
into our code GIZMO, together with state-of-the-art implementa¬ 
tions of MHD in SPH. We systematically compared these meth¬ 
ods to the results from grid codes and analytic solutions on a wide 
range of test problems, and find that the MFM and MFV methods 
are at least competitive with state-of-the-art grid MHD codes, and 
in many cases may have some advantages. 

Critically, we find that our meshless methods can indeed cap¬ 
ture phenomena like the MRI (with the correct growth rates and 
mode structure), formation of magnetically-driven jets in collaps¬ 
ing cores, MHD fluid-mixing instabilities (the Rayleigh-Taylor 
and Kelvin-Helmholtz instabilities), and both sub and super-sonic 
MHD turbulence. In fact we find convergence on these problems 
in our MFM/MFV methods is comparable to, and in some cases 
even faster than, AMR codes using constrained transport. This sup¬ 
ports the similar conclusions found by Gaburov & Nitad orij ( |201 1) , 
studying the MFV method on a smaller set of problems|' |Histor- 


20 To compare performance of the MHD algorithms, consider a simple test 
which essentially counts cycles per second. We initialize a 3D box with a 
standing linear wave of negligible fractional amplitude (10 6 ; just so the 
values of gradients, and solution of the Riemann problem, are not trivially 
vanishing), and evolve it for a short time. We consider a pure hydro case, 
and then add a trace magnetic field that does not alter the dynamics. 


21 We note that our major extension of the work in |Gaburov & Nitadori[ 
4201 1) is to extend the methods to include (and test) the MFM method, as 
well as to conduct a systematic comparison with SPH and grid codes, on a 
wide variety of tests not considered in that paper. We have also made many 
subtle improvements of the MFV algorithm (all described here and in Paper 
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ically, these problems have been difficult for SPH; our new meth¬ 
ods, however, do not suffer from the low-order errors that typically 
cause problems for SPH MHD. But we also show that even SPH, 
with the most current implementation of MHD, is able to capture 
most of these phenomena, albeit at the cost of larger kernels and 
some “by hand” adjustment of artificial dissipation parameters. 

5.1 The Divergence Constraint & Conservation 

Most importantly, we find that, using a state-of-the-art implemen¬ 
tation of the |Dedner et al.| j2002| l divergence-cleaning scheme (re¬ 
discretized appropriately for our new methods), we are able to 
maintain V ■ B « 0 to sufficient accuracy that divergence errors 
do not corrupt the solutions to any of our test problems. Typi¬ 
cally, this amounts to a “worst-case” h |V ■ Bj/|B| ~ 0.01 even 
in highly non-linear problems evolved for long times. In smooth 
flows and/or highly-resolved problems, more typical values are 
/z | V - B|/|B| ~ 10 -4 . 

This is important: without any V ■ B correction, the MHD 
equations are numerically unstable, and most problems will either 
crash or converge to unphysical solutions. The simplest “fix” in the 
literature is to just subtract the unstable terms, the so-called Pow- 
|ell et a 1. | < ] 1 999ft or “8-wave” cleaning. However, we show that the 
Powell cleaning alone converges to the wrong solution on most test 
problems. 

The problem is, in certain types of MHD discontinuities, 
Powell-cleaning alone produces the wrong jump conditions, even 
in the limit of infinite resolution, because the errors occur across a 
single resolution element and are zeroth-order. This has been shown 
before for a limited range of problems in fixed-grid codes; here we 
show the same applies to a wide range of problems in all the meth¬ 
ods considered here. In the Brio-Wu and Toth shocktubes, the shock 
jump conditions are wrong, the same problem leads to qualitatively 
incorrect features appearing in MHD blastwaves (less dramatic ver¬ 
sions of these errors appear in both the Orszag-Tang vortex and 
MHD rotor). In advection of a magnetic field loop, the field strength 
can grow unstably - the same errors disastrously corrupt the non¬ 
linear growth of the MHD Rayleigh-Taylor and Kelvin-Helmholtz 
instabilities, and lead to orders-of-magnitude incorrect growth of 
seed fields in cosmological MHD turbulence (e.g. the Santa Bar¬ 
bara cluster test). In the protostellar jet test, associated momentum 
errors can “kick” the core out of its disk. 

With a good implementation of divergence-cleaning, we find 
that all of these errors are eliminated, provided that converged solu¬ 
tions are considered. This is clearly critical to almost any interest¬ 
ing astrophysical problem. Unfortunately, it means that many previ¬ 
ous MHD studies (see, for example [Dolag & Stasyszyn|2009[[Blir-| 
|zle et al.|2011| |Pakmor & Springel|2012||Zhu et al.|2015[ ), which 

relied only on the simpler Powell-type schemes, may need to be 
revisited. 

It is worth noting that, of the tests we explore here, a combina¬ 
tion of the MRI, protostellar jet launching, magnetic Santa Barbara 
cluster, and non-linear magnetic Rayleigh-Taylor instability, appear 
to be the most challenging to simultaneously capture accurately. We 
encourage authors of future MHD methods papers to include these 
as opposed to only focusing on a subset of problems like the MHD 
rotor, Orszag-Tang vortex, and shocktubes, which we find compar¬ 
atively “easy” and not as useful. 


I); these do not change the qualitative behavior on any tests, but do tend to 
decrease numerical noise. 


5.2 MFM vs. MFV vs. Moving-Mesh Methods 

In all of our tests, we find small differences between our meshless 
finite-volume (MFV) and meshless finite-mass (MFM) methods; 
those (minimal) differences are similar to what we saw in pure hy¬ 
drodynamics tests in Paper I. MFV, with mass fluxes, is able to more 
sharply capture contact discontinuities and minimize overshoot in 
the density/velocity fields around them. However, the additional 
fluxes lead to enhanced “grid noise,” so the method is slightly more 
noisy. 

In practice, the differences are sufficiently small that the “bet¬ 
ter” method will depend on the problem. For some purposes (e.g. 
cosmological simulations), it is extremely useful to maintain ap¬ 
proximately constant particle/cell masses (because the dynamics 
are dominated by gravity in an /V-body solver); this is accomplished 
more elegantly and significantly more accurately with MFM than 
with MFV plus cell splitting/merging (which is more analogous 
to an AMR-type code). But in other cases, high resolution might 
be desired in low-density regions of the flow, in which case MFV 
(possibly used in the mode where cells do not move exactly with 
the fluid velocity) is more natural. 

We have not presented a detailed comparison with moving- 
mesh codes, because a public moving-mesh MHD code capable 
of running the tests here is not available; however, a few of the 
test problems here have been considered in other studies with the 


moving-mesh codes AREPO, TESS, and FVMHD3D 

Duffell & Mac- 

Fadyen 2011; Pakmor et al. 2011; Gaburov et al. 

2012). In each 


of these cases the results are very similar to ours here (especially 
similar to our MFV method). This is consistent with our extensive 
hydrodynamic comparison in Paper I, and expected, given that the 
methods are closely related. 

Most importantly, on all tests both MFM/MFV methods ex¬ 
hibit good convergence properties and capture all of the key quali¬ 
tative phenomena, even at relatively poor resolution. 

5.3 Comparison to SPH MHD 

Historically, it has been very difficult to capture non-trivial MHD 
phenomena with SPH. However, in the last few years there have 
been tremendous improvements to almost every aspect of the basic 
hydrodynamic algorithms in SPH, as well as the specific discretiza¬ 
tion of MHD (see references in § |TJ. As a result, we find that state- 
of-the-art SPMHD is, in fact, able to capture most of the important 
MHD phenomena studied here, including non-linear MRI, dynamo 
effects, magnetic jet launching, and fluid mixing instabilities. 

However, convergence in SPH is still very slow; in almost ev¬ 
ery case, SPH is still significantly more noisy, less accurate, and 
more diffusive at fixed resolution compared to our MFM/MFV 
methods, and requires some “by hand” tweaking of numerical pa¬ 
rameters to give good results on all tests. There are two fundamen¬ 
tal problems: first, SPH requires "artificial diffusion” (viscosity, 
conductivity, resistivity) terms, which are somewhat ad hoc. The re¬ 
sistivity term in particular is challenging, as discussed in § |3.9| the 
correct “signal velocity” and question of whether resistivity should 
be applied at all is much less clear than, say, artificial viscosity, as 
it depends on the type of MHD discontinuity (not just whether one 
is present). We are unable to find a single “switch” the works best 
for all cases, and we show that using an even slightly less-than- 
ideal choice (e.g. using the magnetosonic versus Alfven speeds for 
the resistivity signal velocity) can catastrophically corrupt certain 
problems (such a fluid mixing instabilities and/or jet launching). 
Similar conclusions were reached in |Tricco & Price| l |2013j >. A po¬ 
tential solution to this is the replacement of the artificial resistivity 
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with the full solution of a Riemann problem between particles, as 
in “Godunov SPH” schemes (see |Iwasaki & Inutsuka|201 1| >. Sec¬ 
ond, and more fundamental, SPH has low-order errors which can 
only be suppressed by increasing the neighbor number in the ker¬ 
nel. This leads to an effective loss of resolution and higher diffusiv- 
ity. However not increasing the neighbor number to some very large 
value 100 in 3D) leads to disastrously large errors and noise in 
most of our test problems. Similar problems are well-known in the 
pure-hydro case (see Paper I), but they are much more problem¬ 
atic in MHD, because of how they interact with the artificial re¬ 
sistivity and divergence cleaning terms. Kernel-scale noise seeded 
by the low-order errors produces magnetic divergences, which are 
then subtracted off and damped away, potentially corrupting the 
real solution (and preventing the algorithm from identifying “real” 
divergence errors). For this reason, the typical V ■ B errors and nu¬ 
merical diffusion are ~ 2 orders of magnitude larger in SPH (even 
with >100 neighbors) compared to MFM/MFV methods (with just 

32)0 

Still, provided sufficiently high resolution and large kernel 
neighbor number are used, together with care in choosing the artifi¬ 
cial diffusion parameters specific to the problem, we conclude that 
“modern” SPH MHD can produce accurate solutions. And SPH 
MHD may still have some limited advantages in specific contexts. 
The artificial diffusion operators are wholly operator-split from the 
hydrodynamic operators; when there are extreme energy hierar¬ 
chies between kinetic, magnetic, and thermal energies, it ensures 
that small errors in any of the three terms do not directly appear 
in the others. It can handle free surfaces trivially and maintain nu¬ 
merical stability with vacuum boundaries; however the MHD equa¬ 
tions will not be correct at these boundaries (the zeroth-order errors 
become order-unity, although they are numerically stable). And it 
remains the most computationally simple method we study. 

5.4 Comparison to AMR 

In all cases, our new mesh-free methods (MFM/MFV) appear com¬ 
petitive with state-of-the-art grid-based codes (e.g. third-order PPM 
methods, with constrained transport, and CTU-unsplit integration, 
as in ATHENA). We find no examples where there are qualita¬ 
tive phenomena that either method cannot capture, nor any exam¬ 
ples where we cannot converge to similarly accurate solutions. Of 
course, there are quantitative differences in the convergence rates, 
and errors at fixed resolution, which depend on the method. Despite 
the fact that Eulerian codes can use CT-methods to maintain the di¬ 
vergence constraint, we identify several problem classes in which 
convergence is faster in our Lagrangian methods than in AMR. 

Not surprisingly, these tend to be problems where advec- 
tion, angular momentum conservation, self-gravity and/or follow¬ 
ing large compressions are important - these are the obvious ar¬ 
eas where Lagrangian methods have an advantage. For example, 
we see significantly faster convergence on the field-loop advection 
problem (our MFM/MFV methods produce about the same numer¬ 
ical dissipation as grid methods at 4°-higher resolution, where D 

-- Note that some authors have attempted to control the noise in SPH MHD 
by performing operations only on “re-smoothed” quantities (see |Dolag &| 
|Stasysz yn 2009 Stasyszyn et al. 2013}. This is similar in spirit to increas¬ 
ing the kernel size, and similarly leads to a loss of resolution and increase 
in numerical diffusion. But it is not clear whether the discretized equations 
after re-smoothing actually consistently represent the true hydrodynamic 
equations (they are not, mathematically, the Lagrangian-derived SPH equa¬ 
tions), so it remains unclear whether such methods can actually converge 
(at any resolution) to the correct solution. 


is the number of dimensions.) The mesh-free algorithms are robust 
to arbitrary “boosts,” which degrade the non-moving grid solutions 
on problems with mixing and/or contact discontinuities (e.g. the 
Rayleigh-Taylor and Kelvin-Helmholtz instabilities, Orszag-Tang 
vortex, MHD rotor, and blastwave/explosion problems). The er¬ 
rors caused by these boosts are, of course, resolution-dependent 
(and will converge away in grid codes), but this means that conver¬ 
gence to a desired accuracy on these problems is usually faster in 
our MFM/MFV methods than in grid-based methods, if the fluid 
is being advected at super-sonic velocities. This difference also 
means our new methods are robust to arbitrary velocities in the 
current sheet test (while some stationary-grid methods will crash 
for modestly super-sonic motion around the sheet). Perhaps most 
dramatically, the protostellar core collapse/MHD jet problem com¬ 
bines high-dynamic range collapse, self-gravity, and evolution of a 
global thin disk (angular momentum conservation) - as a result, 
convergence is much faster in MFM/MFV than in AMR meth¬ 
ods. Qualitative phenomena (e.g. the jet momentum/mass) start 
to converge at resolutions as low as 10 4 cells/particles, compared 
to at least 0.3 —lx 10 7 cells in AMR (“effective” resolutions of 
> 20,000 3 )[^]And the resolution demands become more severe in 
AMR if the disk is rotated and/or moving with respect to the coor¬ 
dinate axes: like in Paper I with a simple Keplerian disk problem, 
this requires >> 512 3 resolution in AMR in the disk for good behav¬ 
ior, plus a comparable number of elements along the jet, to prevent 
it from numerically grid-aligning (artificially bending) and being 
destroyed. In contrast, our new methods are trivially invariant to 
such rotations and boosts, at any resolution. 

Of course, on other problems, grid methods converge more 
rapidly. In smooth, pressure-dominated flows, the “grid noise” is 
minimized in truly fixed-grid (non-AMR) methods, so convergence 
in the highly sub-sonic regime (Mach numbers < 0.01) is usually 
faster. On shock-tube type problems (e.g. the Brio-Wu & Toth prob¬ 
lems above), where our errors are dominated by the noise intro¬ 
duced by divergence-cleaning and non-zero V ■ B errors, we see sig¬ 
nificantly faster convergence in grid-based codes that can use con¬ 
strained transport to eliminate these errors entirely (a similar factor 
~ 4° as above). And of course, by virtue of not being Lagrangian, 
in high-dynamic range problems Eulerian codes will better-resolve 
/ow-density regions of the flow. 

In short, we see no “inherently” superior method between 
AMR and our new mesh-free methods, simply differences in the ac¬ 
curacy achievable at fixed resolution or computational cost, which 
are highly problem-dependent. 

5.5 Areas for Improvement & Future Work 

This is a first study, and there are many potential areas for im¬ 
provement. Several possibilities discussed in Paper I (more accurate 
quadrature rules, generalizing to higher-order fluid reconstructions, 
better-optimized kernel functions) apply as well to the MHD case. 

23 It is a common mistake to refer to “kernels” in mesh-free methods as 
“resolution elements” the same way single-cells are referred to in grid 
codes. This is wrong. In our MFM/MFV methods, the correct identifica¬ 
tion is to think of each particle as equivalent to a cell in a grid code (with 
about the same “effective resolution per cell/particle”). The kernel repre¬ 
sents the number of neighbors in causal contact for hydrodynamics; so the 
correct analogy is to the stencil used in a grid code (number of neighbors 
with adjacent faces, plus those needed for gradient calculations). Our de¬ 
fault choice for MFM/MFV, then, of 32 in 3D, is actually quite similar to 
what is obtained in moving Voronoi-mesh, AMR, and higher-order (PPM) 
Cartesian grid codes. 
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The most dramatic improvement to the meshless methods 
here, however, would come from incorporating constrained trans¬ 
port. Recently, jMocz et~aI7[ ( j2014a( l demonstrated that constrained 
transport could be successfully incorporated into moving-mesh al¬ 
gorithms; there is no conceptual reason why the algorithm de¬ 
scribed there cannot be applied to our MFM/MFV methods, since 
they are conservative finite-volume schemes with a well-defined set 
of effective “faces’" and a partition of unity. In contrast, there is no 
clear way to generalize this to SPH (given the inherent zeroth-order 
inconsistency in SPH derivative operators, it is not clear whether it 
is possible under any circumstances to derive a CT-SPH method). 
However, the method in |Mocz et al.| ( |2014a| ) has not yet been ex¬ 
tended to three dimensions and to adaptive timesteps, in a efficient 
manner which can run in competitive time. Therefore we have not 
considered it here, but it is certainly worthy of more detailed explo¬ 
ration in future work. 

Absent a complete CT implementation, some progress might 
be made using locally divergence-free gradient representations, or 
(similarly) vector potentials. Previous efforts have been made in 
this area in both SPH MHD and discontinuous Galerkin methods 
(see e.g. Miyoshi & Kusano||2011[ |Mocz et al.|2014b[ |Stasyszyn| 


|& Elstner||2015 i. These can offer some improvements; however, 

they usually sacrifice consistency and/or conservation, and by only 
providing locally divergence-free terms, it is by no means clear that 
they actually reduce the relevant errors driving numerical instability 
( |Price|20T0l >. But again, further study is needed. 

In SPH, some errors (e.g. the zeroth-order errors) are inherent 
to the method. Others, however, could be decreased. There has been 
considerable work on improved switches for the artificial viscosity; 
similar work is needed for the artificial resistivity (following the 
work of Tricco & Price 2013). In particular, it would greatly expand 
the flexibility of the method if a switch were devised which could 
correctly interpolate between the relevant propagation speeds of the 
resistivity (which depends on the type of MHD discontinuity). 
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APPENDIX A: THE SMOOTHED-PARTICLE MHD 
IMPLEMENTATION IN GIZMO 


As noted in the text, the magnetic terms in our implementation 
of SPMHD follows that in the series of papers by |Tricco & Price] 
< |2012[ [20T3^ . The full SPH hydrodynamics algorithms in GIZMO 
are given explicitly in Paper I. This includes both our implemen¬ 
tation of “traditional” SPH and “modern” PSPH. We note that the 
additions for MHD are independent of whether TSPH or PSPH is 
used, but we will adopt the modern PSPH formulation as our “de¬ 
fault.” 

With the hydrodynamics in place, the additions for MHD in 
SPH are as follows. First, note that we directly evolve the conser¬ 
vative quantities VB and mijj, as in the MFM and MFV methods. 
As noted by |Price & M onaghan| < |2004|) ; |Bate et ah] p014| >, this 
has several advantages over explicitly evolving B and ip - namely, 
it improves the overall conservation properties of the code, elimi¬ 
nates significant errors associated with compression/expansion of 
the particles/fluid elements, reduces noise from poor particle or¬ 
der, allows us to write several of the key equations in manifestly 
anti-symmetric form (allowing conservation to be maintained even 
under individual time-stepping), and greatly simplifies some of the 
hydrodynamic calculations. Overall we find a net improvement in 
accuracy with the VB approach as opposed to directly evolving B in 
SPH; however there can be some advantages to the latter (for exam¬ 
ple, slightly reduced storage and conservation of an initial B x = 0 
in ID MHD; see [Price & Monaghan|2004] >. With this in mind, the 
primitive variables B and i/< are constructed in SPH front the con¬ 
served variables as follows: 


b , s mi s E (VK) , 

Vi mi 

wo, 

mi 

pi = '^m j W (x,- - xj , hi) 


(Al) 

(A2) 

(A3) 


where p, is the usual SPH density estimator (constructed from 
neighbors, hence distinct from the actual density field evaluated at 
the position of particle i in our MFM/MFV methods). 
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The conserved variable (VB) is evolved according to: 
d(V B), 


dt 
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Here x,y = x, — x ( , 12 is the number of spatial dimensions, ® de¬ 
notes the outer product, V ® B is the u x u gradient matrix of 
B, computed using our second-order consistent matrix gradient 
method (and jV®B| is the Frobenius norm of the matrix), x is 
the unit vector x/|x|, v)) st is the fast magnetosonic wave speed be¬ 
tween particles, c Sj , is the particle sound speed (usually computed 
as Cj,,- = ( 7 Pi/pi) 1 ' 2 ), and va is the usual Alfven speed. 

The first term in d(VB)i/dt is the induction equation. The pre¬ 
cise form of this is derived exactly from the SPMHD Lagrangian; 
any other form will introduce errors in conservation and potential 
numerical instabilities (see Tricco & Price 2012). The Q, terms here 
and throughout are derived from the same Lagrangian approach fol¬ 
lowing Springel & Hemquist (20021, and account for variations in 
the “smoothing length” h between particles^] 


The second term in d(VB),/dt is the artificial resistivity (Price 


|& Mo naghan 2005). Just like artificial viscosity and conductivity 
in the pure hydrodynamic case (still present here), artificial dissi¬ 
pation terms are necessary in SPH for all hydrodynamic quantities 
to account for discontinuities. However unlike artificial viscosity, 
artificial resistivity is still needed in rarefactions to prevent numer¬ 
ical instability, so this term is always “active” between neighbors 
(independent of whether they are approaching or receding). The 
form here is motivated by (although significantly different from) 
the dissipation in a Reimann problem; the important aspect is the 
“switch” a B , which one would like to have a large value when there 
is a sharp discontinuity in B, and a vanishing value in smooth flows. 


Tricco & Price 

(2013 

i: a B oc h |V ® B|/|B. In 

Tricco & Price 


<|2013|>, the authors show this is considerably more accurate, and 


24 The functional form of the H terms is slightly different here versus in 
|Tricco & Price|(2012} , because we use the particle number density n,. rather 
than pi, to determine the SPH smoothing length, but this has almost no 
effect on our results in any test problem. The two formulations are identical 
if particle masses are equal, which is also usually the case. 


less diffusive away from shocks, than the “standard” (constant-a s ) 
approachp~| 

Note that, in the artificial resistivity term, the appropriate “sig¬ 
nal velocity,” cf ; , is ambiguous (the physically correct value de¬ 
pends on actually solving the relevant Reimann problem to deter¬ 
mine the type of MHD shock). Tricco & Price ( 2013b adopt the 
mean fast magnetosonic speed, |Price & Monaghan 1 2005} adopt the 
RMS Alfven speed; here we adopt a compromise. When va 3> c s , 
we find the |Tricco & Price| l |2013fr speed usually gives better results 
(the same conclusion they reached in their test problems). However, 
when c s 3> va, and there is particle disorder (either because of mo¬ 
tion induced by external forces or near discontinuities), the problem 
is that the zeroth-order SPH errors always seed non-trivial (~ 1%- 
level) af, so even if there is a smooth, continuous gradient in B, 
the resistivity is triggered and the magnetic fields are damped on 
a sound-crossing time. For some of the problems in this paper, for 
example the MHD RT and KH instabilities and the SB cluster, this 
suppresses the mean field by an order of magnitude, and leads to 
a qualitatively incorrect solution. This is remedied if a wavespeed 
which vanishes with |Bj (e.g. a multiple of the Alfven speed) is 
used. We therefore allow for the us e of either wavespeed in prin¬ 
ciple, with the parameter af in Eq. 


A9 


but adopt af = 1 as our 
“default” (i.e. simply set cf, to the mean Alfven speed). However in 


the tests described in § 3.9 we consider af = 0, i.e. setting cf to 


the mean fast magnetosonic speed. 

The third term in d(VB)i/dt is the divergence-cleaning term, 
ocVV>. The particular functional form is again Lagrangian-derived; 
as pointed out in |Tricco & Price (2012), this is especially impor¬ 
tant, since not just any form of the gradient estimator can be used. 
Rather, it is necessary to use one which operates in appropriate con¬ 
jugate pairs with the gradients used for the V ■ B estimation and 
pressure-gradient (hydrodynamic force) operations, or else the re¬ 
sulting cleaning scheme will be numerically unstable, and simply 
fail to clean the “correct” divergences. 

The conserved variable (mtp) is evolved according to: 


d(mip)i _ 2 
dt ~ Vs ‘ 6 


1 ah (TF E m > ( B ' “ B ')' ^Mij(h/) 

^ M Pi 


( i \ a P Vsi 8> 1 

- f—r 

J kern ”/ 


(All) 


The first term in d(imp)i/dt is the corresponding source term for 
the divergence-cleaning field (the hyperbolic term). Again, the 
functional form inside the summation is strictly tied to the func¬ 
tional form of the cleaning in d(VBW<i?rjThe second term in 
d(rml))i/dt is the parabolic damping rj Here /k e m is a constant de¬ 
fined for convenience that depends on the kernel shape (= 1/2, 1/3 


25 We actually further improve on this formalism, by using our matrix- 
based gradients to determine V ® B. Just like with the higher-order artificial 
viscosity switches proposed in|Cullen & Dehnen (2010|, the use of second- 
order consistent gradients (as opposed to the zeroth-order inconsistent SPH 
gradient estimator) greatly improves the accuracy of the switch (helping it 
trigger in the “correct” locations). 

26 Here we follow the “difference” formulation from Tricco & P rice I 
<20l2> , which they show provides the greatest stability and minimizes er¬ 
rors among the formulations they compare. 

27 Note that, in closer analogy to grid-based methods, we evolve ( tmp ), and 
not, for example, (Vip) or -ip. This produces an essentially identical set of 
equations for the evolution as in |Tricco & Price|{2012) . We do not include 
their additional advection term 

- E'_ T. m i ,n i ( v ; -D) ■ ViWij(hi) (A12) 

2 UiPi 
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for the cubic/quartic splines), for h defined as the kernel radius of 
compact support. 

In both of these, v S i g ,i is the maximum signal velocity calcu¬ 
lated between all neighbors, as described in the text. This signal ve¬ 
locity is modified compared to the standard hydrodynamic case by 
the replacement of the sound speed c S] , with the fast magnetosonic 
speed vj) st between particles. That replacement applies for all places 
where the signal velocity appears - for example, in the artificial vis¬ 
cosity terms, and calculation of the CFL condition/timesteps. 

The only remaining equation is the magnetic force: 


d(m\)i\ 


clt 


=Yj m 

i 

—B, y nti nij 


Vi yi. 

h ■ ViWij(hi) + ■ V,W u (hj) 


n.pj 


n iPj 


H/Pi ^jPJ 


M, = 


B,-®B,- ]B,- 

Mo 2 P‘0 


I 


(A 13) 
(A14) 


where M is the Maxwell stress tensor and I is the identity ma¬ 
trix. The first term here is the usual MHD acceleration d\/dt oc 
p~ l V ■ M; again, the particular functional form of the gradient op¬ 
erator derives necessarily from the SPH Lagrangian (see Price &| 
Monaghan 2005; note it is essentially identical to the form of the 
Lagrangian-derived SPH hydrodynamic force in Springel & Hem-| 
|quist|2002] with M replacing the pressure P). The second term is 
the Bprve et al.|( |200l| implementation of Powell 8-wave cleaning 
- namely, subtracting the unphysical part of the equation of mo¬ 
tion proportional to V ■ B. As in the main text, this is necessary to 
prevent catastrophic numerical instability (here in the form of the 
tensile instability). 

There are now four numerical parameters that must be set: the 
artificial resistivity terms a® in , Qmax> and the divergence-cleaning 
terms o p , a /,. The divergence-cleaning parameters are discussed ex¬ 
tensively in the text; we find a best compromise on all problems in 
this paper using the “default” values oi, = 1, a p = 0.1 PI As noted 
there, cr/, is a “nuisance” parameter that can be folded into the defi¬ 
nition of the cleaning speed - we include it here only for complete¬ 
ness. For the artificial resistivity terms, unless otherwise specified 
we take a® in = 0.005, a® ax = 0.1. We have experimented exten¬ 
sively with these, and find these are best compromise values. The 


in the evolution equation for ip in |Tricco & Price] (20 1 2f . Like them, we 
found that this term does nothing to improve behavior on our tests, and 
can de-stabilize the cleaning procedure in simulations where the velocity 
divergence is large (e.g. cosmological runs), without additional timestep re¬ 
strictions. However we have run almost every test in this paper with the 
term active and find (provided proper care is used in timestepping) that the 
differences are negligible. Also following |Tricco & Price | (20 12^ , we have 
experimented with an artificial dissipation term for ip. However, because 
in SPMHD there is no ip flux, we find (as these authors did) that this pro¬ 
duces no improvement in performance on any test problems here, and only 
increases the numerical diffusion. 

28 While o p ~ 0.1 is typical for mesh-based codes and appears optimal 
for our MFM/MFV methods, |Tricco & Price||20T2| favor oy, = 0.1 in 2D 
but o p = 0.8 — 1 (i.e. more rapid damping of ip) in 3D for SPH. However 
this is based just two tests; and in at least one o p ~ 0.1 actually minimizes 
the maximum value of h | V • B|/|B|. Moreover different definitions of c/, 
make direct comparison difficult. While we certainly agree that larger o p 
is beneficial on some tests, we find it can lead to substantially larger diver¬ 
gence errors on others; hence we adopt the more “conservative” cleaning 
parameter ( o p ~ 0.1). 


a max = 0-1 choice follows Tricco & Price (2013 ; a much larger 
value (e.g. a® ax ~ 1) completely diffuses away the fields in several 
of our test problems (e.g. the RT instability, MHD rotor, SB clus¬ 
ter, protostellar core collapse) and dramatically over-smooths shock 
jumps in others (e.g. the Zeldovich and Toth problems), leading to 
systematically incorrect solutions. But if a^ 0.1, the method 
is incapable of properly capturing strong, magnetically-dominated 
shocks. The lower limit is less important; a® in <0.1 is important 
to prevent excess diffusion that suppresses field growth in e.g. the 
proto-stellar disk problem, but a® in > 0 greatly reduces the noise 
and post-shock oscillations that seed low-order SPH errors. 

This is sufficient for SPMHD. However, in running a large 
suite of shocktube tests (at the suggestion of the referee), it be¬ 
came clear that our default implementation of artificial viscosity 
(described in detail in Paper I and taken from the “inviscid SPH” 
prescription in |Cullen & Dehnen|2010|l is not ideal for some MHD 
problems. Most dramatically, in the Brio & Wu (19881 shocktube, 
with small neighbor number (“SPH-lo”), adopting the prescription 
from Cullen & Dehnen (20101 (or the slightly modified forms in 
|Hu et al.||2014 or Paper I) leads to very large oscillations in the 
post-shock velocity over the domain where the internal energy is 
large. As shown in Paper I, on pure-hydro problems the method be¬ 
haves well. The issue in MHD appears to be insufficient viscosity in 
regions with sub-sonic noise when the accelerations are primarily 
transverse. The simplest solution is to enforce a constant artificial 
viscosity, but this seriously degrades the performance of SPH on 
other problems (hence the reason for these switches). After some 
experimentation we adopt the following compromise between ex¬ 
cessive diffusion and noise: the functional form of the artificial 
viscosity follows the hydro case in Paper I (Eq. FI6), except we 
replace the sound speed with the fast magnetosonic speed and in¬ 
crease the minimum viscosity from a m in = 0.02 to a m m = 0.05. As 
well, the dimensionless parameter aois set by 


atmp — ' 


0 ((<f[V - x]/dt)i > 0 , or (V-v), >0) 


C^max | (r/[V ■ \]/dt ), 
|(d[V ■ y]/dt)i\ + (c/h) 2 


(otherwise) 


a 0 ,,(f + At) = < 


Olrnp (atmp A O!0,,(f)) 

a tmp + (ao, ,(t) - a tmp ) A ' 

(atmp < ao,/(f)) 


(A 15) 


where in Paper I we took c = 0.7 c s ,,- and h = fk C m h,. Here, we use 
h- -A C/ken, hi + (A 16) 

c —*■ 0.2 (c~l + ||w|| -1 ^ _1 (A 17) 

where ||v|| denotes the Frobenius norm of v and ||v|| = ||v||. The 
prescription from Paper I has the effect that the viscosity van¬ 
ishes very quickly whenever the mean velocity gradient is resolved 
and/or the compressive accelerations are sub-sonic. With the mod¬ 
ifications above, a velocity field which has large fractional noise or 
rate-of-change of the velocity divergence (relative to the velocity 
itselD will also “trigger” viscosity (even if the flow is sub-sonic or 
some of the noise is transverse). This is sufficient to significantly re¬ 
duce the noise in the |Brio & Wu|(T988j ) problem, while having only 
small effects on almost every other problem here. It does, however. 
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somewhat degrade performance (via larger numerical viscosities) 
on some problems like the rotating Keplerian disk in Paper I. More¬ 
over, as written here, this term violates Galilean invariance; we find 
this can be restored with similar (slightly more noisy) results by in¬ 
stead using || v|| = a v MAX;(||v,- — v,-||) where MAX; refers to the 
maximum among neighbors and a v ~ 10. 

This completes the SPMHD implementation. 


APPENDIX B: ON THE USE OF VARIABLE 
WAVESPEEDS IN DIVERGENCE-CLEANING 
OPERATORS 


In |Dedner et al.|{2002| >, the authors showed that their divergence¬ 
cleaning method is numerically stable and guaranteed to reduce | V ■ 
B|. However, strictly speaking, their proofs are valid only if the 
wave and damping speeds c* and r (r = c 2 p /c\ in their notation) 
are constant in both time and space. 

In that case, if we neglect any other fluid forces and evolve the 
system only under the influence of the mixing terms in B and ip, 
the evolution equations for ip and V ■ B take the form: 


d 2 ip 1 dip 
~W + r~di 


£v 2 ip 


d 2 (V ■ B) 1 d(V-B) 
dt 2 + r dt 


-c?V 2 (V- B) 


(Bl) 

(B2) 


i.e. both ip and (V ■ B) obey a damped wave equation. 

However, as discussed in the text, it is desirable to vary Ch 
and r. Imposing a constant Ch equal to some global maximum 
wavespeed essentially forces all elements (cells/particles) to use a 
constant, global timestep (a tremendous numerical cost), because 
they must all satisfy the Courant (CFL) condition for the wave 
equations above (otherwise | V - B| can grow unstablylp^Moreover 
in many problems the “fastest” wave speed may reside in regions 
which have wildly different physical properties and are causally 
disconnected from others. This also leads to wavespeeds which dra¬ 
matically exceed local physical signal velocities. 

If we allow Ch and r to depend on space and time, following 
the same derivation as in|Dedner et al.|(|2002|l gives: 


d 2 iP + 1 


dt 2 t(x, t) dt 


9 ^ -cl(x,t)'V 2 ip = 


(B3) 


dt 


dip ip\ 1 dc 2 h , d 
+ F) cl^-^dtx r 


9 2 (V-B) 1 a(V-B) _ 2/ _ 


dt 2 


- + 


t(x, t) dt 


-cl(x,t)V-(V-B)= (B4) 


ipV 2 ^J+ 2 ViP-V^j + 

(V ■ B) V 2 (cl) + 2 [V(V ■ B)] -V(cl) 


Now the equations have time-and-space dependent coefficients, and 
source terms dependent on derivatives of c/, and r. And unfortu¬ 
nately, these extra terms cannot be eliminated by_simply modifying 
the original source terms in the 


Dedner et al. 


1 2002) scheme □ 


29 The CFL condition for Eq. |b 1 |B2| requires all elements have timesteps 
A t{ <C r ~ MIN(/iy)/MAX(v s i g j) (the criterion in c/, is slightly less de¬ 
manding). This is equal to or smaller than the normal minimum timestep 
for any single element in the si mulation, MIN(Af, ) ^ MIN(/z 7 /v s i gj7 ). 

30 Consider e.g. modifying the |Dedner et al.| (2002) assumption that the 
correction term scales as dB/dt = —\7if with the insertion of an arbitrary 
function / such that dB/dt = — /(x, t,if, (V • B)) 'Wif, or adopting alter¬ 
native operators T>, g defined such that T>(if) +g(V • B) = 0 (in their for- 


So consider the extra terms. First, take Ch = Ch(t),r = r(t), i.e. 
the coefficients are spatially constant, but change in time. This cor¬ 
responds to the original implementation proposed by |Dedner et al.| 
( |2002| ), and most subsequent work. This is generally not a problem. 
The time-derivative terms appear only in the if equation (Eq. |B3[ ), 
which has no direct physical consequence. Eq. |B4| for V • B remains 
a damped-wave equation, but with time-dependent coefficients. So 
long as their time variation is sufficiently slow, one can apply the 
usual Wentzel-Kramers-Brillouin (WKB) approximation and show 
that V • B still behaves as a damped wave. Because Ch is chosen 
to be a maximum wavespeed in the domain, and l/r some max¬ 
imum damping rate ~ c/ 7 /MIN (/*,), the variation of Ch and r in 
time (which evolve with the physical properties of the system such 
as p, B, etc.) will always be slow compared to the local evolution 
timescale for the damping wave (provided the system obeys the 
CFL condition in the first place), and this is easily satisfied| 3 4 

Similarly, if Ch and r depend on position, we will recover the 
desired behavior so long as c/,(x) and r(x) are sufficiently smooth. 
From Eq. |B4| the corrections terms do not change the behavior of 
the system if | Vr|/r and | V c/ ? |/c/ 7 are < | V(V • B)|/| V • B|p| 

We expect and show in the text that (since it is a local numer¬ 
ical error term) | V ■ B|/| V(V • B)| ~ h\\ so our choice of c/ 7 and r 
(which depend on local physical properties) should guarantee this 
condition is satisfied wherever the flow is resolved^However, cau- 


mulation, g = 1, T> = cf/ 2 d/dt + r —1 ). If we do this and attempt to recover 
Eq. |B i| one can show that the correction terms can simply be folded into 
a new field if' that obeys Eq. |B2| with constant coefficients. If one desires 
(V • B) to obey a damped wave equation with constant coefficients, and to 
build the “correction terms” out of linear operators acting on (V • B) and 
any arbitrary field if, then the |Dedner et al.|(2002) scheme is the most gen¬ 
eral possible solution. 

31 We caution that if the time-variations of c/ 7 and r are large on a wave¬ 
crossing time (e.g. c~// ] \dch/dt\ > Ch/\L\, where |L| ~ |V • B|/|V(V • 
B)| hi), the WKB approximation is invalid and |V • B| can converge 
to a constant or grow. If for example c/ 7 evolves on some dynamical time 
?dyn? the requirement for stable behavior is c/, > hj/t^ n . This is another 
reason to choose c/, to be the maximum (local) wavespeed (which always 
satisfies this). Choosing a slower wavespeed, even if uniform in space, can 
de-stabilize the cleaning. 

32 More precisely, a ID analysis following |Dedner et al.|(2002) gives the 
following sufficient (although not strictly necessary) criteria for stability of 
the divergence-damping: 


dx(c 2 T) 


< 1*1 , 




L /l 

dx(F h r) 


<l Y 




, _ a-(v-B) 
(V-B) 


\k\ 2 C 2 T 2 


1*1 


2 | k\ 2 c 2 h r 2 


„ _ d 

d x = — 

dx 


(k < 0) 

1*1 (* > 0 ) 


(B5) 

(B6) 

(B7) 


i.e., the sign of d x [c^T (V • B)] (and d x [c^ (V • B)]) must match the sign of 
<9. v (V • B). This is particularly easy to see if we consider just the parabolic 
term in a ID case (where any deviation from B x = Bq = constant represents 
V • B f 0). In this case we obtain the heat equation: d t B x = d x (c^r d x B x ), 
and deviations from Bq are diffused away, provided the condition above 
is met. Physically, the stability condition corresponds to the requirement 
that the spatial dependence of c/, and r does not introduce new local ex¬ 
trema in c|rV • B that are not present in V • B. These conditions are ob¬ 
viously satisfied (in 3D) if we have |Vc||/c^ and \'Vcj i r\/cj j T much less 
than |V(V ■ B)|/|V - B|, our (more restrictive) criterion above. 

33 Note that c/, appears in two places: in the Riemann problem (the first- 
step update of the normal component of B) and in the source term for if. 
In the former, we are solving the one-dimensional Riemann problem in 
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tion is needed when the flow is poorly-resolved, especially in multi¬ 
phase, turbulent systems (where physical properties vary rapidly). 
By linking c*,, and t; to the maximum signal velocity among in¬ 
teracting neighbors, as opposed to speeds at i alone, we maintain 
smoothness even with kernel-scale noise. 

As a result, we confirm in our tests that | V ■ B| does not grow 
unstably to levels that swamp the physical solutions. However, one 
could probably improve the cleaning (especially in noisy, multi¬ 
phase flows) by doing more to ensure c/,(x) and r(x) remain suf¬ 
ficiently smooth. For example, we have experimented with calcu¬ 
lating Ch,i and Ti following our standard method in a first particle 
sweep, then computing a kernel-weighted average (c*},-, (r), on 
a second sweep; this ensures smoothness on a super-kernel scale. 
Alternatively, we have considered an effective slope-limiter which 
limits the magnitude of c 2 h l and 77 such that among interacting 
neighbors, the particles which are local (kernel-scale) extrema in 
V ■ B remain so in cl V ■ B and cl r V ■ B. Although there are hints 
of minor improvement on a couple of problems, we do not find 
these make a large difference and have not experimented with them 
extensively. However, the issue merits investigation in future work. 


Of course, much too large a value compared to other characteristic 
speeds in the problem will mean that p> cannot grow (and therefore 
cannot remove divergences), while much too small a value can lead 
to a “lag” in the response of the cleaning to V ■ B. 

We have therefore considered a variety choices: 


1 MAX 
c r,i — ^ Dig,; 


c t, i — /kem hi T" fastest 


TfasLt = MAX, 


2 /kem hi / k , 


~~r~ (o,; + V X,) 

xn hi ' ' 


1/2 


C t ,;=MAX 


1 MAX 

2 V sig ,i > v /Vm 5 ^fastest 
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APPENDIX C: COSMOLOGICAL INTEGRATION OL 
DIVERGENCE-CLEANING 

The modifications in our code necessary for cosmological integra¬ 
tions are described in detail in Paper I. These are, for the most part, 
unchanged. As described therein, it is easy to show that if we define 
appropriate co-moving units, the cosmological expansion is auto¬ 
matically handled, and the Riemann problem is locally unchanged, 
provided we convert into physical variables before solving it. This 
is identical in MHD. For the magnetic field B, we consider the 
co-moving field B r = a 2 B, where a = 1/(1 + z) is the usual scale 
factor. This is invariant under pure adiabatic expansion in a Hub¬ 
ble flow in ideal MHD. Since we evolve the conserved variable 
(VB), and the volume/length units are also co-moving, we simply 
have (VB) C = a^ 1 (VB). For the divergence-cleaning field ip. the 
“correct” co-moving units are slightly more ambiguous; how ip be¬ 
haves in a Hubble-flow expansion actually depends on the ratio of 
timescales (whether, for example c s > va, because ip depends on 
the fastest wavespeed). We have experimented with simply assum¬ 
ing ip c = a 3 ip (the appropriate choice in the typical cosmological 
case, when either c s or the inter-particle fluid velocity |v| is much 
larger than va), or ip c = a 5 ^ 1 if’ (appropriate if va dominates), or ex¬ 
plicitly solving for the evolution terms from expansion; in practice 
we find this makes no detectable difference to any problem, since 
the physical t/>-field growth and decay time to respond to evolving 
V ■ B (~ h/ch , where h is a resolution element and c* the fastest 
wave speed) is always vastly shorter than the cosmological expan¬ 
sion/Hubble time. 

APPENDIX D: THE DAMPING SPEED 

In § [2] we note that the divergence-wave is damped with a source 
term where usually 77 = hi/(a p c T .i) with c Tj , the damp¬ 

ing wave speed. There is no a priori obvious choice for this speed, 
and the ability of the scheme to damp divergence does not depend 
sensitively on the choice, provided it satisfies the conditions in §[B] 

operator-split fashion, so the only states that matter are the two faces - we 
can either enforce a constant C/, within each Riemann problem (our default 
approach described in the text) or explicitly solve for the two-c;, system (de¬ 
scribed in §[F] since the waves propagate away from the discontinuity at the 
face into media with constant c/, on either “side,” there is no instability). 


The first is the standard signal velocity v^f/2 (this is the de¬ 
fault choice in SPH MHD). We find this works well in every prob¬ 
lem here, although in some cases it produces excess dissipation of 
the fields because it operates too slowly for some subset of parti¬ 
cles when ip is very large (since the speed at which ip can induce 
changes in B is not accounted for with this velocity). 

The second (v/^,i) is the particle-based fastest-possible mag- 
netosonic speed, with the additional term 2 ipi/v^f. This term can 
be thought of as representing the “potential magnetic energy” in the 
ip field - while usually negligible compared to the Alfven speed va, 
it is certainly possible, on some problems, that ip grows until it 
reaches values |t/>,| |B|/|v S i g | - this is clearly in the limit where 
damping of the ip field is desired. Therefore we find this choice 
works much better than the pure magnetosonic speed. This choice 
works comparably well to v^f, but is less than ideal in some cases 
where the particles have super-sonic local approach velocities (not 
accounted for here). 

The third choice (Vfastest) is similar to the default choice in 
|Dedner et al.| ( [2002j >. namely the fastest wavespeed and/or signal 
velocity in the entire domain. This works well on idealized test 
problems - including almost all of the tests in this paper; how¬ 
ever, it makes little sense for high-dynamic range problems like 
cosmological or galaxy/star formation simulations. In those cases 
the medium is highly multi-phase, so there is a huge range of local 
fastest wavespeeds, often in regions which are not even in causal 
contact. We therefore find that this produces too-efficient damp¬ 
ing (hence less-efficient cleaning) in the slowly-evolving regions 
of these problems. 

The fourth choice (fastest) is similar to the maximum 
wavespeed Vfastest, but instead sets 17 directly to the minimum across 
all particles (independent of whether they have small/large local 
volumes hi). In a uniform-grid code this is identical to Vfastest- Here 
we find it works comparably well, again in idealized test problems, 
but has the same problems in inherently multi-scale problems. 

The final choice, therefore, represents our best attempt at a 
“compromise” between these. We take the maximum of either the 
signal velocity v^f/2, the local magnetosonic speed (plus ip) 
Vf^^, and some multiple with e* <C 1 of Vfastest- We adopt this as our 
default for all problems here, with ei, = 0.01. However we stress 
that all of our qualitative results are robust to any of the choices 
above for c Tp - we find only minor quantitative differences (in 
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many of the test problems here, these are completely indistinguish¬ 
able). 


APPENDIX E: ADDITIONAL FLUX-LIMITERS FOR 
DIVERGENCE-CLEANING OPERATIONS 


In Paper I we describe our slope-limiting procedure for reconstruc¬ 
tion; the same is used here. However, in a couple of cases (e.g. the 
isolated galaxy disk and Santa Barbara cluster), some additional 
flux limiters greatly help in improving numerical stability. These 
are specific to the divergence-cleaning (ip) terms, so do not directly 
affect our reconstruction of the physical quantities. They will alter 
how efficiently the V ■ B errors are cleaned; while this is important, 
it is not the only consideration for stability. 

In particular, numerical instability can arise owing to the 
“mixing” terms between ip and B which appear in the updated val¬ 
ues B' x i j and ipij in the Riemann problem (see Eq.|l4|. In the limit 
where there is large particle disorder (e.g. poorly-resolved, turbu¬ 
lent shocks) and the particles are being rapidly re-arranged by non- 
MHD forces, or if boundary particles are only able to find a couple 
of neighbors (if vacuum boundaries, which are not recommended 
for MHD, are used) then one can have \ip\ 3> i’^? g AX |B,j and the im¬ 
plicit instantaneous update of B (~ (ipL — iPr)/ci,) in the Riemann 
problem can be unstable (since small residuals in ip which are not 
completely damped by the time particles locally re-arrange can lead 
to large changes in B); this is related to the discussion in § [b] 

One solution would be to simply drop these terms. However, 
this prevents the divergence-cleaning from acting across single¬ 
particle discontinuities, which (as discussed in the text) can lead 
to incorrect jumps. We find a more accurate, robust, and flexible 
solution, which works well for all problems in this paper, is to sim¬ 
ply apply an additional set of limiters for ip. 

In the Riemann problem, we modify Eq.[l4]to be: 


K. ,J = \ ( s ', l + b' x R ) + (ip L - ip R ) 


2 Cft i 


'Ipij = 7z('lpL + 1pR)-\-'lpij 


ipfj = ^ [B' X) l B' x R ) 


a^'ij = MIN 


, 0 Ch,ij \B X L +B X s | 

1, OL^jj 


\IpL-IpR 


(El) 

(E2) 

(E3) 


This is a flux-limiter on the implicit ID Riemann problem between 
B' and ip at a discontinuity (or equivalently, we can think of it 
as limiting the slope of the ip discontinuity). The coefficient ct/p 
should be < 1 for stability; the precise value is (like all limiters) 
set by a balance between stability and diffusion. Our experiments 
prefer a!J, = 0.75. For the source terms for B which are oc (VX7ip)* 
in Eq. ED we also modify 


(VVVO,* X' ' - A - -»• (YV1>)lo + al,i(VVip)* B 
(Wip)U = - 

(VViP)*B = -J2^ Ai J 


y— (tpL ~\~ 1pR)jj A 

- b AlJ 


a B ^i = MIN 1, 

c _ j d(VB), 
dto 


IOC 2 

\{VViP)* B \\ 

2 | i 0.5v sig ,,-CEB), | 2 


(E4) 

(E5) 

(E6) 

(E7) 

(E8) 


where d(VB)i/dta represents the value of d(VB)j/dt calculated 


for element i including all other fluxes and source terms ex¬ 
cept [Wip)/ B . Similarly, in the source term for ip (Eq. Ill, 
dip/dt oc (W • B)f clu we limit the effective value of (VV ■ 
B)* allowed to a maximum = 100 V) |B,|//z,-. The pre-factor 
is of course arbitrary but should be > 1. Finally, we check 
each timestep whether | ip/\ > av max |B,j with a = 10 2> 1 and 
Vmax = MAX(0.5 v sig ,i, Vfastest, J c] t .); if it exceeds this value, 

we impose dip/dt = MAX(0 , dip/dt) (ip < 0) or dip/dt = 
MIN(0, dip/dt) (ip > 0). This just corresponds to increasing the 
(already arbitrary) ip- damping rate super-linearly when ip becomes 
very large. 

Even in the Santa Barbara and galaxy disk problems, these 
limiters almost never act. Usually, when they do, the particles are 
in a situation (e.g. at vacuum boundaries) where the fluxes are un¬ 
resolved and should not, in any case, be trusted. Therefore, the fact 
that this limiting procedure allows somewhat higher V ■ B errors 
(by making the i/>-based cleaning less aggressive) is a small price 
to pay for maintaining numerical stability. 


APPENDIX F: A TWO-WAVE FORMULATION OF THE 
DIVERGENCE CLEANING TERMS IN THE RIEMANN 
PROBLEM 

As discussed in the text, in Eq. [14] the normal component of B 
and the divergence-cleaning term ip are implicitly updated accord¬ 
ing to the solution of an independent one-dimensional Riemann 
problem, before solving the MHD Riemann problem. The solution 
given there assumes a single wavespeed Ch.ij = MAX [vf,L, Vf.«] 
(the maximum of the fast magnetosonic speed on left and right 
sides of the problem) for the divergence-cleaning wave at the dis¬ 
continuity between the left and right states. Following a solution 
proposed by E. Gaburov (private communication), we could instead 
assume two independent wavespeeds, cl — Vf,r and c R = Vf.R on ei¬ 
ther side of the discontinuity. This yields the solution 

B x ij = -■- [ ClB x r + CrB x r + ipR — ipR^ (FI) 

c L + c R l 

ipij = -■- [c* IpL + c L ipR + ClCr (B' x L -B' x R )] (F2) 

cl + Cr 

This trivially reduces to the solution in Eq.[l4]for cl = Cr. We have 
considered this formulation, instead of the default in the text, for all 
test problems in this paper. In all cases, the differences are small. 
The two-wave formulation introduces some additional dissipation 
and/or grid noise, depending on the problem, however it is also 
more stable in situations with large particle disorder (essentially 
because it provides an up-wind weighting of the ip and normal-B 
terms). 
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